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CHAPTER I 


INTRODUCTION 

Ground wind induced flow fields around surf ace; < . ^ i 

obstructions such as buildings, bridges . and; o their man-made i 
structures have long been of , interest in structiairal design.;; :; 
For many years wind loading has been; :the subject of studies 
in various wind tunnels. Most of these studiSiS, however, 
have been related to specific problems rather than to the;. . 
more general and systematic flow field investigation which 
‘can. .provide a code of practice for the building engineer. :.i 

More recently such fundamental investigations have 
attained neW; importance through environmental considerations; , ; 
as well as through the many aerodynamical problems encoun- , 
tered in the design of airports for V/STOL aircraft or 
helicopter ports in large metropolitan areas. Operating low 
speed aircraft between buildings in regions of steep velocity 
gradients or large fluctuations through vortex fields or ■ 

recirculation zones is very hazardous. A clear under- 
standing and detailed knowledge of the flow field around 
typical buildings is therefore a necessary source of 
information to minimize the dangers of these problems. The 
need for analytical methods to predict local atmospheric 
motions influenced by buildings or similar bluff surface 
obstructions is then obvious. 



In formulating an analytical model to completely 
describe this complicated flow situation, the complete 
equations of motion must be considered. However, practical 
methods for carrying out a solution of such equations are 
limited to: nimerical approaches which are presently in the 
development stage. , . : ; ; : 

Some of the more recently developed methods which 
were quite successful shall be discussed here in brief. 
However, only very few methods incorporate the effects of 
turbulence. Most of them are applicable to laminar time- 
dependent .two- or three-dimensional incompressible viscous 
flow problems with various boundary conditions. One of 
these methods, which uses pressure and velocities as primary 
variaibles , and which can be applied to confined; flows as 
well as flows having free surfaces is the "marker-and-cell” 
(MAC) method developed by Harlow and Welch [1, 2]^ at Los 
Alamos . Certain marker particles are introduced into the 
flow calculation to indicate the respective fluid configu- 
ration. The flow equations which are described in Eulerain 
space (full Navier-Stokes equations, incompressible); are 
formulated in finite difference form, in both space and time 
variables. The derivation of these finite-difference 
equations is based upon the following sequence of events by 

^Numbers in brackets refer to similarly numbered^ . 
references in the Bibliography. 
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which the whole configuration is advanced from one step to 
the next . 

1. The velocity distribution is known at the 
beginning of the cycle either from an initial 
condition or from the previous cycle:. 

2. The corresponding field of pressure is> Ji,-; 
calculated . 

3. The components of acceleration are calculated? 
from these changes in velocity are computed .and 
added to the previous values. 

:.j: ; 4. The marker particles are moved to new posit i<ons 

according to the velocity determined for their 
locations . 

The method is capable of satisfying a free-slip or 
no-slip condition at a rigid wall and a pressure boundary 
condition at a free surface. Computations have indicated 
considerable numerical stability and comparisons with some 
experimental results showed very close agreement. However , 
finite difference approximations introduce truncation errors 
that can obscure the effects of real viscosity and influence 
the stability of the solution. If only the mean effects of 
turbulence on a flow are of interest, such difficulties can 
usually be avoided, because the effective turbulent vis- 
cosity is often larger than the molecular viscosity and is, 
therefore, not likely to be obscured by finite difference 
errors. A method of this type, proposing ZIP differencing 
[3] is presented by C. W. Hirt [4] . However, ZIP 
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differentiating cannot be readily incorporated into the 
marker-and-cell technique, since its variables are defined 
at different mesh points. 

R. S. Hotchkiss [5] quite successfully employs an 
extension of the MAC method for his calculations of three- 
dimensional flows of air and particles around structures. 

The full nonlinear Navier-Stokes equations for incompres- 
sible flow with buoyancy are solved. The effect of turbu- 
lence is resolved through a constant eddy viscosity. The 
method is capable of satisfying five different kinds of 
boundary conditions, e.g., inflow, outflow and periodic ‘ *-■ 

boundaries, along with rigid walls that have either free- 
slip or no-slip conditions. In addition to the Navier- 
Stokes equations, the time dependent heat equation is solved 
to incorporate features of thermal buoyancy (Boussinesq 
approximation) . The computer program also includes the 
simultaneous solution of a transport equation in order to 
determine distributions of particulate matter. The results 
are displayed as three-dimensional perspective views drawn 
by a computer. Similar equations, however, using a variable 
eddy viscosity have been used by J. W. Deardorff [6] for his 
numerical investigation of the idealized planetary boundary 
layer. The nonlinear equations of motion are solved for the 
ideal case of neutral stability, horizontal homogeneity of 
all dependent mean variables except the mean pressure. No 
buoyancy effects are considered. A discussion and details 
of the numerical method are given in [7] . ' 
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Although it was possible in all these examples to 
study full three-dimensional problems with methods such as 
the MAC method, these calculations are still extremely time 
consuming and require the largest computers now available. 

A simplified marker-and-cell; .method , (SMAC) .has been given), by . 
A. A. Amsden and F,. ,H. Harlow .[8]. Another extension ;Of the: < 
SMAC method is proposed by P. I. Nakayama and N . , C i , ;Romero 
[9] for the solution of incompressible transient flows that 
are almost three-dimensional. The equations of motion are 
two-dimensional but contain functions that account for the 
effects of a third dimension. The method may also be used 
for the simple incorporation of internal obstacles in two- 
dimensional flows. , 

Besides the MAG method and its numerous modifi- 
cations, some; other numerical experiments with a difference 
model for the time-dependent Navier-Stokes equations have 
been reported by Schoenauer [10]. He concludes that the r 
space mesh size must be inversely proportional to the 
Reynolds number. A coarse net produces a numerical "turbu- 
lence" which tends to the physical turbulence as the mesh 
size goes to zero. However, the existing computers are not 
fast enough for fluid flows which have Reynolds numbers in 
the turbulent regime. 

Similar to Hotchkiss [5] , D. Djuric and J. C. Thomas 
[11], apply a. numerical model (following Harlow and Welch [1] 
and .l?earjdorff [6] ) to atmospheric boundary layer calculations, 
especially to transport and diffusion of gaseous air 
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pollutants in the vicinity of tall buildings. The equations 
used in this model follow those derived by Ogura and 
Phillips [12J in which the "anelastic" approximation is 
utilized. In these equations the basic density p is 
constant in the horizontal planes and varies only in height. 
The Coriolis term is not included because of the small scale 
of motion under consideration. Periodic boundary conditions 
were used except for pressure, which had a slight decrease 
in wind direction to provide a driving force which moves the 
air . 

Ro B. Lantz and K. H. Coats [13] use a quite 

different model for their three-dimensional calculation of 
spread and dilution of air pollutants. Their mathematical 
models are similar to the "Gaussian plume" models [14, 15, 
16] and their extensions which can include topographic 
effects [17, 18] . A numerical solution of the three- 
dimensional material balances for pollutant flow and for the 
air stream is given. The pollutant material balance 
requires the solution of the three-dimensional convection- 
diffusion equation. Wind flow over uneven surfaces or 
around simple structures is calculated by numerically 
solving an equation for the velocity-potential which is 
modified such that it includes a variation of the horizontal 
velocity with height (logarithmic or power law) . The 
numerical schemes used in the solution of the two balance 
equations are said to be efficient ones requiring minimal 
computer time. The mathematical description of individual 
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eddies downstream from obstacles is felt not essential to 
the pollutant dilution problem. The increased dispersion of 
the pollutant due to eddy formation is approximated by 
increasing the eddy diffusity downstream of the obstacles. 

There has long been a difference of opinion as to 
which of the possible forms of the equation of motion is the 
most suitable for niimerical solutions. Some workers, as for 
example, Harlow and Welch [1] , prefer to retain the veloc- 
ities and pressure as dependent variables. Others, such as 
Aziz and Heliums [19] feel that it is advantageous to use 
instead the vorticity and stream function as dependent 
variables. Gosman, et al. [20] , share the latter opinion. 
Their finite-difference model for steady two-dimensional 
flows, as described by elliptical partial differential 
equations, has proved to be very successful in a large 
number of problems to which it was applied. The model is 
capable of handling turbulent mean flows with variable fluid 
transport properties. Several turbulence models such as the 
Prandtl mixing length concept or a turbulence kinetic energy 
model using the Prandtl-Kolmogorov hypothesis have been 
employed. The computational procedure consists of solving 
the differential equations for the fundamental conservation 
laws of mass, momenttim and energy which are complemented by 
auxiliary relations for the transport properties. Because 
of the established versatility of the model and because of 
the generality of its framework, the model shall be used in 
this investigation to study the turbulent flow with an 
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atmospheric velocity profile over a bluff surface obstruc- 
tion. 

An alternative approach, approximating atmospheric 
motibns over surface obstructions has been carried out by 
Frost, et al. (1973) [21], in extending the concepts of 

boundary layer theory. Their method was applied to the 
specific geometry of a semi-elliptical cylinder. The 
characteristics of atmospheric shear flow over a rough 
terrain (disturbed boundary layer) are coupled with the 
turbulent boundary layer equations using Prandtl's mixing 
length hypotheses. Two approaches are presented to inco:?-'’^ 
porate the pressure field and boundary conditions which 
exist within the large viscous region over the obstruction. 
The first considers a region in the immediate vicinity of 
the body in which the pressure distribution and outer 
boundary condition on the velocity are computed from poten- 
tial theory for flow over the elliptical cylinder. The 
second approach considers a much larger region of influence, 
extending from the surface to the undisturbed flow at large 
heights above the obstruction and uses a vertical pressure 
decay function to blend into the undisturbed flow region. 

The main conclusions drawn from the study were: 

1. Localized maxima in wind speed occur at the top 
of the semi-elliptical obstruction. 

2. An increase in the elliptical aspect ratio 
decreases the wind speed within the^ boundary ' 
layer at the top of the ellipse and .returns it 
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to the logarithmic distribution characteristic 
of the undisturbed flow. 

3. Increase in surface roughness affect the flow by 
decreasing the velocity in the boundary layer, , 
with the most pronounced effect occurring near . 
the surface of the smaller aspect ratio ellipse. 

4. Reynolds number has a negligible effect on the 
overall flow for the range of Reynolds numbers 
considered. 

To investigate the validity of the Prandtl mixing 
length theory for atmospheric flow in disturbed regions 
Frost and Harper [22] extended the above approach, still 
solving the boundary layer equations with equivalent boundary 
conditions, but using the conservation equation for turbu- 
lent kinetic energy (TKE) instead of Prandtl ' s mixing 
length (PML) to model the effect of turbulence. Comparing 
the respective flow fields predicted by the two models of 
PML and TKE it is found that there is only a small differ- 
ence in the predicted velocity profiles, primarily in the 
regions of strong pressure gradients, caused through the 
diffusion and the convection of turbulence kinetic energy 
not accounted for in the PML concept. However, while the 
PML model does not provide any information about the turbu- 
lence structure of the flow field, the TKE model gives 
physically quite meaningful values of turbulence intensity 
levels which are presented and discussed in the study. 
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Although the method described has been used in both 
cases for the flow over semi-elliptical geometries only, it 
can generally be applied to any two-dimensional body shape. 

In the course of the following investigation it shall there- 
fore be used for the calculation of the flow field over a 
bluff surface obstruction. Again the Prandtl mixing length 
concept will be employed in view of its reasonable results 
and its numerical efficiency. The potential model, however, 
used to compute the pressure gradient term and outer 
boundary condition for the flow over the semi -elliptical '• 
geometry, will be replaced by a more sophisticated model Of>'" 
inviscid shear flow over a bluff body. 

This inviscid model, which in itself already repre- 
sents an interesting analytical solution to the problem, is 
described in Chapter II. Its governing equations are 
developed and the required empirical input is discussed. 

The application of this model to atmospheric flow conditions 
is outlined and some typical results are presented. 

Chapter III deals with the aforementioned boundary 
layer approach incorporating the pressure field input from 
the inviscid model of Chapter II. The governing turbulent 
boundary layer equations, their numerical solution procedure 
and the initial-and boundary-conditions are briefly reviewed. 
Solutions are given for selected parameters characterizing 
the undisturbed atmospheric velocity profile. 

In Chapter IV the problem is then approached by 
solving the Navier-Stokes equations following the method of 
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Gosman, et al. [20] . Three different turbulence models are 
presented and compared — the PML concept, Prandtl-Kolmogorov' 
TKE model with a given turbulence length scale distribution 
and a TKE model with a simultaneous solution of a transport 
equation fp.r the turbulence length scale. In contrast to ■ 
the boundary-layer approach where the equations of motion 
are solyesd in terms of the primitive variables u, w and p, 
the governing equations are written in terms of the vor- 
ticity u and the streamfunction \p and the other physical 
properties involved. The equations are then expressed in a 
commcgii'i form, differing only by a source term peculiar to the 
property which the equation represents. The transformation 
of the differential equations into difference equations 
which are solved by successive substitution, the incorpo- 
ration of turbulence into the model as well as the 
respective boundary conditions used are explained. Finally, 
some results for different parameters of the approaching 
wind profile are presented and discussed in the reniainder of 
the chapter. , 


\ 



CHAPTER II 


INVISCID FLOW MODEL 
I. INTRODUCTORY REMARKS 

A suitable model for flow over a bluff body should 
at least include the occurrence of flow separation and a 
wake formation. For inviscid flow models this leads auto- 
matically to a free-streamline approach. However, the 
classical assxamption of the free-streamline analysis (23J„ is 
that of a quiescent fluid on the inner side of the free,j.^. 
streamline such that the free streamline remains one of 
constant pressure and velocity. This results in a constant 
pressure wake of infinite extent which is an inadequate 
representation of the actual case. The model yields a base 
pressure Pj^ which is equal to that at infinity, whereas it 
should really be Pj^ << p^, thus leading to negative base 
pressure coefficients Cpj^. The conventional free-streamline 
theory therefore has to be applied with moderation. 

A somewhat different free-streamline theory for an 
inviscid shear flow over a bluff body, featuring an upstream 
separation bubble and a closed downstream wake region with 
variable pressure, shall now be described. 

Parkinson and Jandali [24] developed a simple theory 
for a two-dimensional incompressible uniform flow external 
to a bluff body and its wake. The desired flovj separation 
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points are made the critical points of a conformal trans- 
formation to a region in which a surface source creates 
stagnation conditions at the critical points. The stagna- 
tion streamlines then transfoann to tangential separation 
streamlines in the physical plane (Figure 2.1). The 
position and the strength of the source is determined by the 
empirical requirements of the separation position and its 
pressure coefficient. 

This model was then used by Kiya and Arie [25] and 
extended for their calculation to a separate flow past a 
bluff body attached to a wall on which the approaching 
turbulent boundary layer has been replaced by a hypothetical 
inviscid uniform shear flow (Figure 2.2). In addition, it 
utilizes a solution found by Fraenkel [26] to incorporate 
the formation of a corner eddy in front of the body. This 
model admits analytical solutions and automatically yields 
closed streamlines in front of the body which are geometri- 
cally very similar to those observed in practice. Like the 
Parkinson-Jandali model it includes a finite wake width and 
a pressure distribution on the separation streamline which 
decreases asymptotically towards the corresponding free 
stream value at infinity. The theory requires input of four 
empirical parameters which depend on the geometry of the 
body. These parameters are: 

1. Location of the front face stagnation point. 

2. Stagnation point pressure. 

3. Location of the separation point. 



Figure 2.1. Physical- and basic transformation-planes with separation points 
Si and Sj as critical points [24]. 
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4. Separation point pressure. 

In the present analysis the Kiya and Arie [25] model 
has been extended further. The original model yields a wake 
region which is finite in width, however, infinite in length. 
Through the addition of a sink equal in strength to the 
existing source, at a suitable downstream distance in the 
transfojrmation plane it is possible to close the wake and 
thus eliminate the unrealistic feature of the original rear 
model region. This involves an additional empirical 
parcuneter, i.e. , the location of the rear reattachment point. 

II. DEVELOPMENT OP GOVERNING EQUATIONS FOR 
INVISCID FLOW OVER A FENCE 


The details of the inviscid flow model shall now be 
developed for the case of the flow over a fence. 

We consider a two-dimensional incompressible, 
inviscid, steady uniform shear flow past a bluff body ASB as 
shown in Figure 2.3. The velocity profile at a large dis- 
tance from the body shall be 

u* = U* + k*y* (2.1) 

where U* is the dimensional velocity at the wall and k* the 
dimensional vorticity 


k* 


9u* _ 3v* 
3y* “ 9x* 


( 2 . 2 ) 


which is constant throughout the flow field. Introducing a 
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stream function rp* by 


u* 


d\b* 

dy*' 


V* 


ax* 


( 2 . 3 ) 


where u* and v* are the velocity components of the velocity 
vector q* in the x* and y* direction respectively/ i. e._ / 


q* = u*i + v*5 


( 2 . 4 ) 


The continuity equation 
div q* = 0 

a terfw 

is then automatically satisfied. The vorticity vector J2* is 
given by 


= curl 5*,- 


( 2 . 5 ) 


and from Equation 2 . 2 


V^ip* = k* 


( 2 . 6 ) 


We can nondimensionalize the coordinates by a reference 
length h* 


X = 


X* 

h*‘ 


= Yl 

h* 


( 2 . 7 ) 


and the velocities by a reference value U*gf 


u = 


u* 


U* 

ref 


V = 


V’ 




Uo = 


fi*' 

_ Up 








( 2 . 8 ) 


•1 B.: 
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and thus rewrite Equations 2.1 and 2.2 as 

Uq^ = Uo + ky 

and 

_ 9u 3v 
^ "^y “ 3x 

Furthermore , Equation 2 . 6 becomes 
= k 
where 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


” o - * - !*- h* <2-12) 

ref 

is the dimensionless stream function. It can easily be 
shown that Euler's equation of motion for steady flow 

2 

V ^ - qx curl q ® -grad p (2.13) 

or 


“It-'' I? '-IS 

“M-“W=-|S <2-24) 


simply reduces to Equation 2.11 which is therefore the 
governing equation of our flow problem. 

On the boundaries, we require the following condi- 
tions to be satisfied: 
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\|» * const. - (2.15) 

oh the wall and solid surface and 

a 

Uo + ky; H 0 (2.16) 

at large distances away from the body. If we subdivide the 
stream function into two parts, i.e., 

j ky^ + y 

then 

» k + 

and from comparison with Equation 2.11 we get 

» 0 (2.19) 

as the new governing equation with the transformed boundary 
conditions 

y + j ky^ » const. (2.20) 


(.2.17) 


(2.18) 


on the wall and solid surface and 


3y 


U 


0 » 


9x 


0 


( 2 . 21 ) 


at large distances away from the body. 

Through the introduction of a function $ which is 
related to by the Cauchy-Riemann relations 
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3$ d'V 3$ 3Y 

J 

we can define a complex potential 
W{z) = «• + iy 


( 2 . 22 ) 


'ly. 


(2.23) 


where 
z = X + iy 



(2.24) 


describes our physical z-plane. We can now transform the 
wall boundary A^A and together with the fence contour 
ASB into the upper half of a new plane called the ^ -plane 
(see Figure 2.3, page 17) where ' " 


^ ^ + XT\ ■ ::'^5 r.Mv ; - - ^ (2.25) 

and i ■ 


- ..+5^ h 2a 

by using a suitable transformation 
z'= f(c) 


(2i26) 

.O ' J o 


I 


(2.27) 


Under the assiamption that this transformation, behaves as 
z ~ kiC + kjJlnc + ks + 0(?“M (2.28) 

P .i . -an. i % ao r.^puc.: ••■i -.u 

Fraenkel [26] has found' a solution to Equation 2.19 which 
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satisfies the boundary conditions for Equations 2.20 and 
2.2i/ r.e./ 

'¥ = Im(Wj + Wjj) (2.29) 

where 


Wj = Uo • ki • C 


(2.30) 


and 


W 


II 



z!ii4Ldr 

? - c 


(2.31) 

f -i. f 


In order to have a free streamline which originates from the 
separation point S, extends in the downstream direction and 
reattaches at R, we have to add to our stream function V a 
combination of sinks and sources which are appropriately 
located. For this purpose a third complex plane, the Z- 
plane is introduced by the transformation 

? = Z + ^ (2.32) 

Figure 2.3, page 17, shows this plane with its 
source-sink arrangement. In addition to the source at 
X = ae^ employed in Reference [25] a sink is placed down- 
stream at the location X = ae . 

The complex potential of the source-sipk system in 


this plane is given by 
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W_tt = — [An(Z - ae^) + An(Z - ae”^) - iln(Z - ae*^) 
III 7T 


- Jln(Z - ae”*^) + C] 


(2.33) 


If we add this to Equation 2.29 we get the streaun function ip 
of the resulting flow as 


t|> = i ky^ + Im(Wj + + Wjjj) (2.34) 

The complex velocity 

w=u-iv (2.35) 

;i ) 

can then be determined as 


u - iv = ky + ^ + i ^ (2.36) 

or with the Cauchy-Riemann relations 


u 


iv = ky 


+ 


dW 

dz 


which is 


u - iv = ky + 


dWj 


dW 


II 


dW 


■ar 


III 


dZ 



(2.37) 


(2.38) 


Since the angle of intersection of the curve at S is not 
preserved in the z-plane the point S is the critical point 
of the transformation, i.e.. 


dz 


f ' (0) = 0 


(2.39) 
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In order that the velocity at this point remains finite/ we 
require from Equation 2.38 that 



dW, 


II 


dW. 


Ill 


d? 


dZ 


dZ 


?=o 


0 


(2.40) 


The three complex planes are related by 


-j ^ V. 




( 2 . 41 ) 


, a :. ■: ( 


and 



f 





i it : 


J i > . 



To use Praenkel's solution (Equation 2.19) one has to 

.n. ; :■= 

establish a relation between the 2 - and C-plane from 

Equations 2.41 and 2.42 and find the constants ki and k 2 in 
Equation 2.28. Kiya and Arie [25] have found and to 

be^ ' - 


.. ..n ' n : 

’) i:.' ■ -1 0 :•'=:? H V UC' \ , ;! 

and 


. (2.43.),, 

U. : 


“ii - rirH) <2.44) 

However, to close the separation region, had to be 

determined differently as already given by Equation 2.33. 

.fi' < I ;5' pn. ■ r: i f- ;■) i 

One can divide it into three parts 



r 


'^iii ”iiip ^ 

"hi = - Pj(Z)l + C 

where 

Pg(Z) = An(Z - ae^) + £n(Z - ae"®) - £n Z 
F^(Z) = £n(Z - ae*) + in(Z - ae”^) - An Z 


25 

(2.45) 

(2.46) 


(2.47) 

(2.48) 


In texins of these three complex potentials, the complex 
*ve^pcity in the physical plane can he expressed by Equation 

>. i 

2.38, where from Equations 2.41 and 2.42 


dZ _ Z^ 
dc ” Z2 - a^ 


and 


dz _ Z^ + 
dC “ Z2 - a2 


(2.49) 


(2.50) 


Note, the latter expression has a Simple zero at the sepa- 
ration point Z = ia. Q is now deteirmlned from Equation 2.40 
as 


Q = 2ira 




cosh 3 * cosh 6 
cosh 6 - cosh 8 


(2.51) 


Substituting the foregoing results into Equation 
2.34 and rearranging source terms gives the following 
expression for the stream function 


I 
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ij, = 


1 

2 


ky* + Im 





+ |fFg(Z) “ FgCZ)] 


+ C 


(2.52) 


The complex velocity from Equation 2.38 


u - iv “ ky + 


■ar ”35“ * ~3T~ 


dZ 

dCj 


• il 


becomes 


u - iv = ky + 


(u. + l|a) + I . c . 


C - 2a 
C + 2a^ 


+ Q. z 

TT Z ^ - a ^ 


f 2(Z 

[z^ 


- a cosh 3) 


2aZ cosh 3 + a^ 
2(Z - a cosh 5) 


Z 2 - 2aZ cosh 3 + a2 


Z^ - a^ 
Z”^' " a'2 


(2.53) 


?,A 


• "'i 


A'- 


(2.54) 


Since the nondimens ional vorticity k is constant throughout 
the flow field/ Euler Equation 2.13 can be integrated to 
give Bernoullis equation 



-kip — 


const. 


(2.55) 


or 


+ p . ^ 


(2.56) 


allowing the determination of the pressure coefficient 



r 


Cp = - (u^ + v^) - 2k (Ij)^ - Ip) 

which reduces on the body surface to 

p — ' 

^Psurf ° surf i 
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(2.57) 


(2.58) 


Assvuning that separation occurs at a given base pressure 

1; wn:: v; j . 

Cpj^, assumed constant over the whole rear side of the fence, 
we can write ' s j’' 

Cpb = - v| (2.59) 

As is required through Equation 2.40, that the separation 

velocity be finite:, Vg is deterraihfe^i Equation, 2.54 by 

^ - 

iv_ =-lim(u - iv) 

^ ; .. . i .. 

Z ^ ia ■ ^ , • 1 ' ; ; 

y •> iZa 

' ' • . ' i ' > . : " I y 

C ^ o ‘ ' ; ' ( 2 o 60 ) 


By LiHopital's rule; 


_ ( ' 4kal cosh 3 + cosh ~6‘ 

® ' cosh 3 • cosh 6 , ; 


III. REQUIRED EMPIRICAL INPUT 


(2.61) 


An inviscid shear flow model for the flow over a 
fence, yielding a realistic upstream separation bubble and a 

finite wake region has been formulated in the preceding 

i , ’I ^ ■ 

section. The input of five empirical pareuneters xs required 

to complete the model. These are: 


I 
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1. The front face stagnation point location 

z = iv 
max ■'max 

2. The front face stagnation point pressure 
^Pmax 

3. The separation point location 

4. The separation point pressure 
^PS " ^Pb 

5. The dovmstream reattachment point location 

“ ^R 

Using Equation 2.54 relationships between parameters 1, 3 
and 5 respectively and Equation 2.58 for the remaining two 
parameters, one obtains five equations to solve for the five 
unknowns 

Uo , k, 8 , 6 and Vg 

IV. DISCUSSION OF APPLICABILITY OF FENCE 
MODEL TO RECTANGULAR BODIES 

The flow over a rectangular block (see Figure 2.4) 
can be treated similarly to the one over a fence, since the 
external flow field is qualitatively the same as that for a 
fence provided reattachment does not occur on the roof. The 
upstream separation bubble is only insignificantly different 
and the downstream separation must occur at the sharp 
leading edge and thus at the same location as in the case of 
the fence. A different empirical input for the separation 
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base pressure and the downstream reattachment point would 
therefore account sufficiently for the new flow situation. 
Some typical base pressure coefficients for different 
building geometries are shown in Table 2.1. It is pointed 
out that the rectangular block buildings generally have 
higher Cpj^ values than that currently used for the case of 
the fence. 


V. RELATION TO ATMOSPHERIC PLOWS 


Since the inviscid flow model is dependent on 
empirical inputs, it is reasonable to postulate that if 
these inputs are characteristic of bluff geometries in the 
atmospheric boundary layer, then the resulting solutions 
should reflect flow characteristics of the atmosphere. 

One of the empirical input parameters discussed 
previously, which seems amenable to this approach, is the 
base pressure coefficient Cpj^. According to experiments by 
Good and Joubert [27] , pressures on the upstream face of a 
normal plate located in a smooth-wall-bbundary layer are 
determined by a wall similarity law of the form: 


'Pb 


u; 

2 /> 

u* ^ 

ref 



152 log 


h*u^ 


V’ 


- 147 + P(j) 


h* 

6 * 


'J 


(2.62) 


Here P is constant for turbulent boundary layers with zero 
pressure gradient and <j> is a universal function of h*/6* 
which is tabulated in their paper. When h*/6* is less than 



TABLE 2.1 


F 


TYPICAL BASE PRESSURE COEFFICIENTS FOR 
DIFFERENT BUILDING GEOMETRIES 

^ ^ ^Pb 


1.0 0 

1.0 0.4 0.4 

1.0 1.0 1.0 

1.0 4.0 4.0 


-0.85 

- 0 . 6 * 

-0.5* 

-0.3* 


*Source: Peter Sachs# Wind Forces in Engineering 

(New Yorks Pergamon Press, 1972). 
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0.5, as is the case for most atmospheric boundary layers, 
4»(h*/6*) is negligibly small. Thus, 


^Pb “ " 


’ ' 

2f 

[“Jefj 

• 


h*uj 

152 log -rrs 147 


V’ 


(2.63) 


Now Cp^ obeys a wall similarity law in the sense 
2 

that Cpjj (U* ^^/u J ) can be described as a function of h*uj/v* 
only, for the case of the smooth wall. 

Applying this to the case of the rough wall, where 


-vI-* if 

the similarity law becomes: 


(2.64) 



152 log 


h* 


- 147 


(2.65) 


Figure 2.5 shows a graphical representation of this equation 
for several dimensionless surface roughnesses. Extending 
the equation still further in order to relate to atmospheric 
flows and introducing a geostrophic drag coefficient 

Cg » qT (2.66) 

with G* being the geos tropic wind. Equation 2.65 can be 


written as: 



Base Pressure Coefficient 
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^Pb 




(2.67) 


From Csanady [28] , the geos trophic drag coefficient C is a 

f * • z ^ ^ 

function of the Rossby number Bo- « — i.e. , 


= f(Ro) 


(2.69) 


and hence one can rewrite Equation 2.65 as 


'Pb 


G* 




f (Ro) 152 


log ^ - 


147 


(2.69) 


It is proposed that the inviscid solution can thus 
be related to the prevailing atmospheric condition by using 
a Cpj^ value as an input, which is predicted from Equation 
2.65 with known values of uj and z* or from Equation 2.69 
with known values of G, Ro and z*. 


VI. RESULTS OF IKTVISCID MODEL 


The aforementioned model was applied to the fence 
problem using the following empirical input-parameters 

1. Front face stagnation point location 

2. Front face stagnation point pressure 

Cp-ax * 0-5 (2-70W 

3. Separation point pressure 


^PS “ ^Pb 


(2.70c) 
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4. Downstream reattachment point location 

Xj^ = 13.0 (2.70d) 

Conditions 1 and 4 were taken from Good and Joubert [27] . 
Condition 3 was determined from Equation 2.65 or Figure 2.5, 
page 33, for the values 


u 






ref 


0.0625 



(2o 70e) 


These four empirical input parameters ^Pmax' ^Pb 

Xj^ lead to the following four calculated parameters: 

Uo = 0.7071 
k = 1.0245 
cosh 6 = 1.3132 

cosh 6 = 10.725 (2.71) 

The streamline pattern computed on the basis of these values 
is shown in Figure 2 . 6 and an enlargement of the upstream 
separation region in Figure 2.7. As there are only limited 
experimental results available, especially from full scale 
studies, the quality of this flow field description cannot 
readily be determined. However, when compared with the small 
scale experimental results of [27] very good agreement is 
found, as can be seen in Figure 2.8. 

Velocity profiles of the u-component for various x- 
stations are shown in Figure 2.9. Here the simplification 
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of the model to a constant vorticity flow field becomes 
apparent in the almost linear velocity profiles. Neverthe- 
less, the region of retarded flow upstream, the recircu- 
lating separation bubble in front of the fence and the 
accelerated flow region in the vicinity of the fence are 
qualitatively correct. This is also the case for the plot 
of the velocity gradient along various streamlines in 
Figure 2.10. 

In turn, the computed pressure distribution is 
reasonable. The almost exact agreement of the pressure 
variation along the front face of the fence with experiments, 
which was already found for the infinite wake model [25] , is 
maintained for the modified closed wake flow model. The 
pressure variation along given streamlines near the ground 
and along the separation streamline is shown in Figure 2.11. 
This pressure disl^ribution will be used as the imposed 
pressure field in the viscous turbulent boundary layer 
approach to be described subsequently. 

Relating the model to atmospheric flow conditions by 

varying the base pressure coefficient Cpj^ according to 

relation 2.65 but holding the remaining empirical input 

parameters constant, the results displayed in Figure 2.12 

are achieved. It is found that the maximum height h^,„ of 

max 

the recirculation region behind the obstacle decreases with 
growing friction velocity u* but increases for larger sur- 
face roughness Zq. 






Friction Velocity 



Surface Rouglmcaa <• 


Figure 2.12. Maximum height of recirculation region in dependence of the 
parameters of the approaching wind. 
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Figure 2.13 shows the predicted velocity increase 
due to the presence of the obstruction calculated at an x- 
station of six obstacle heights downstream. represents 

the velocity profile which would exist at this particular 
location in the absence of the obstruction, i.e., the 
approach velocity profile as defined by Equation 2.1. The 
profiles for the velocity increase are given for various 
base pressure coefficients, which can be related to the 
atmospheric conditions with Equations 2.65 through 2.69. 

In conclusion it can be said that this rotational 
inviscid model for the flow over a fence, featuring an up- 
stream corner eddy and a finite wake region is in good 
agreement with experimental results in the vicinity of the 
obstruction. It describes quite realistically the size of 
the front separation bubble, predicts the pressure distri- 
bution along the front side of the fence and determines the 
downstream separation streamline enclosing the recirculating 
wake region. Through the simplification to a constant 
vorticity flow field, however, the far field representation 
which returns to the assumed linear velocity as y -> <» is, 
with the exception of the streamline pattern, not realistic. 
One expects a logarithmic velocity profile in the atmosphere 
and hence both the upstream velocity profile and the velocity 
profile at y " should have this logarithmic form. To 
achieve more realistic far field velocity profiles a 
boundary layer approach which will give a better approxi- 
mation to the flow in a turbulent natural atmosphere is 
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developed in the following chapter. The results of the 
preceding analysis r however, become extremely useful in 
defining the imposed pressure distribution required in 
boundary layer analysis. 



CHAPTER III 


ANALYSIS OF ATMOSPHERIC FLOW OVER A BLUFF 
SURFACE OBSTRUCTION BY THE TURBULENT 
BOUNDARY LAYER APPROACH 


To apply the concepts of boundary layer theory to 
two-dimensional atmospheric flow, consideration is given to 
a homogeneous terrain of infinite extent on which a bluff 
surface obstruction or fence of height h* is located. The 
inviscid flow field around such a fence was established In 
the previous chapter. Far upstream of the obstruction the 
viscous turbulent atmospheric motion is described by the 
logaritlimic velocity distribution 


u* ® 





(3.1) 


where uj is the friction velocity defined in terms of the 
surface shear stress t* 



and ic the von Karman constant of value 0.35 to 0.4 


I . GOVERNING EQUATIONS 

The governing turbulent mean-flow boundary layer 
equations for steady flow within the atmospheric boundary 


47 



48 

l^yer as derived In [21] are: 


au* . 3w* 


0 


(3.3) 


u* 


3u* . 


w* 



(v* + £*) 


3u* l 

32*J 


(3.4) 



(3.5) 


where overbar denotes an ensemble average and is a 
constant density consistent with an adiabatic reference 
State. The pressure p* represents the difference between 
the total pressure and the hydrostatic pressure. It should 
be recalled that, in addition to the conventional boundary 
layer approximations, two other assumptions underlie the 
above equations: 

1. The atmosphere is neutrally stable. 

2. Coriolis effects are negligible, which is a 
reasonable assumption for the atmospheric 
boundary layer below 30 to 50 m, see for example 
Tverskoi [29] . 

The eddy viscosity, e*, may be related to the mean 
flow through the Prandtl mixing length hypothesis in the 
following way: 


e* 


£*2 


du* 

dz* 


(3.6) 


with the mixing length £* as 



Jt* = (zj + z*)k 
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(3.7)’ 


It was already pointed out in Reference [21] that this 
assimption is not strictly correct in regions of large flow 
curvature, but it should suffice in view of further assump- : 
tions imposed later in this report on the curvature of the 
flow and the coordinate system used in the investigation. . 

II. INCORPORATION OP THE INVISCID SOLUTION 

When analyzing conventional boundary layer flow, the ' 
pressure gradient in Equation 3.4 is approximated by the 
pressure variation along the zero streamline determined from 
the inviscid flow solution for the respective body. This 
approximation is justified through Equation 3.5. Typical 
pressure changes for the inviscid flow over a fence are 
shown in Figure 2.11, page 42. When the pressure gradient 
along the zero streamline is introduced into the boundary 
layer equation for the bluff body, however, it may cause 
flow separation which cannot be handled with the boundary 
layer approximation. By successively introducing the pres- 
sure gradients of streamlines further away from the body, 
one can find the first streamline for which the corresponding 
pressure gradient does not yield an upstream separation. 

This streamline will be called the nonseparating streamline 
Ip . The pressure gradient along this streamline is assumed 
to drive the flow in the present problem. This is believed 
to be a reasonable assumption in view of the fact that the 
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primary interest of this investigation is directed toward 
the solution farther away from the obstruction. 

Another input from the inviscid solution into con- 
ventional boundary layer analysis is that a free-stream 
velocity at the outer edge of the boundary layer is pre- 
scribed. In the present approach however, the internal 
boundary layer produced by the surface obstruction has to 
merge with the undisturbed atmospheric boundary layer and 
its logarithmic velocity profile at a sufficiently high 
altitude. Letting the pressure gradient decay to zero in 
the vertical direction merges the velocity profile smoothly 
with the logarithmic one. Therefore, the form of the pres- 
sure distribution introduced into the boundary layer 
equations in previous approaches [21, 22] was given by 


dp* (x* , z* ) _ f dp* (x*) 1 
dx* [ dx* 

^ns 


(3.8) 


where q(z*) is the vertical decay function. An initial 
approximation of this function was given by the following 
second order quadratic decay function: 


q(2*) « 




for < 0.5 


for > 0.5 


(3.9) 


In the current investigation this somewhat arbitrary 
approximation has been replaced by a pressure decay derived 
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from the invlscid model described in Chapter II. It must be 
noted, however, that the vertical pressure distribution 
based on a linear velocity profile will not behave the same 
in the far field as a logarithmic profile. Therefore, it is 
assumed that the respective vertical decay functions are 
proportional to their momentum flux, i.e.. 


q(z*>iog _ 

U| (z^ " U? „ (z* 

log lin 


With the decay of the linear solution given as 


(3.10) 


C^(z*) - Cp 

* Cp(*„3) - c” 


(3.11) 


the resulting pressure distribution finally becomes: 

_ ... .. .2 

dp* (x* , z* ) 
dx* 


dp* (x* ) ' 

‘ ''Pcc 

^log 

dx* J 

, ■ crop ) ■- c~ 

iPns P ns Poo 

^lin 

V . 


(3.12) 


III, COORDINATE SYSTEM 

Since the pressure force driving the boundary layer 
flow is determined along the inviscid streamlines over the 
fence, the coordinate system must also be oriented along 
these streamlines, resulting in the orthogonal system shown 
in Figure 3.1a. The curvature of this coordinate system is 
small throughout the flow regime, except in the vicinity of 
the upstream stagnation point and the downstream reattachment 



(b) Assumed Numerical Coordinates 


Figure 3.1. Coordinate system used in boundary layer 
approach for flow over a fence. 
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zonCf where the slope of the zero streamline is dis- 
continuous. It is argued, however, that at a location 
somewhat removed from the region of large curvature the flow 
will be driven by the pressure distribution along the non- 
separating streamline. This is in agreement with Hunt [30J 
as described later. Alternately one may look at the 
approximation as being an analysis of flow over a solid body 
defined by the nonseparating streamline which envelops the 
upstream and downstream separation regions. Since the 
primary interest in the present report is the flow field 
above the obstruction where aircraft operations occur, the 
boundary layer approximation as posed here is expected to 
provide meaningful results. Since the boundary layer con- 
cept is generally a first order approximation to viscous 
flow, neglecting the higher order effects produced by 
curvature, it is not expected that the present boundary 
layer analysis will provide accurate results in regions of 
strong curvature, as for example, in the upstream vicinity 
of the obstacle. On the other hand, there is little gain in 
transforming the boundary layer equations into curvi-linear 
coordinates thereby introducing additional complications. 

As a consequence the calculations have been carried out in 
the Cartesian coordinate system (x*,z*) shown in Figure 3.1b. 
The x*-direction is measured along the inviscid streamline 
with the z*-axis extending perpendicular to it at each x*- 
station. The resulting velocity profiles calculated in this 
coordinate system are then assumed to exist perpendicular to 
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the zero streamline in the physical plane. This is believed 
to be a reasonable assumption for the following two reasons: 
first, the good geometrical agreement of the separation 
streamline of the inviscid solution with experimental 
results, as can be seen from Figure 2.8, page 38, and, 
second, that these geometrical effects of the recirculation 
region on the flow now enter the equations through the 
imposed pressure gradient. 

Hunt [30] has shown that a first order analysis can 
be justified if 


H*_ Sin{6*/z* 
L* &n(L*/z^) 


(3.13) 


where H* and L* are the characteristic height and length of 
the disturbance and the thickness of an internal boundary 
layer calculated from 

6 * 6 * 

lF ° zy “ (3.14) 

In the present approach H* and L* would be chosen as shown 
in Figure 3.2, with H* being the maxlmiom height of the 
separation region h* as given in Figure 2.12, page 43, and 

lUclX 

L* being the horizontal length of the combined upstream and 
downstream separation region, also specified by the inviscid 
solution. For a typical surface roughness of z* = 0.005 h* 
one would then get from the above condition Equation 3.13 


with 3.14: 
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0.1693 « 0c6497 
or for 2 j = 0c05 h* 
0.1693 « 0,5858 


which are reasonably satisfied. 

The no-slip condition implied in the boundary layer 
solution at the lower boundary causes the calculated 
velocities to vanish along the zero streamline. This 
becomes physically unrealistic where the zero-streamline 
separates from the surface to become a free or dividing 
streamline o Therefore, the following model for calculating 
velocity profiles through the recirculating wake region 
behind the fence is proposed in form of an approximate 
solution using an integral technique under investigation by 
Kaul and Frost [31], Figure 3,3 illustrates the flow 
regions considered and how they are matched with the present 
boundary layer profile. 

A polynomial expression 

u*(x*,z*) « A(x*) •+ B(x*)z* + C(x*)z*^ + D(x*)z*® (3.15a) 


is assumed with the following governing boundary conditions: 


1. u*(x*,0) * 0 


at z* = z* 


(3.15b) 


2 . 


a^u* 


1 = 1 ap* 


3z*2 V*p* X 


at 2 * = 


2*=0 


2 * 


(3.15c) 


numerical Boundary Layer Solution 
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3. 



9z* 


at z* = z^ 


(3.15d) 


4. u*(x*,z*) 



at z* = z* 


(3.15e) 


5. 



u*(z*)dz* = 0 


(3.15f) 


Thus there are five equations for the five unknowns Zj, A, 

B, C and D. Condition 1 implies the no-slip condition at 
the wall. Condition 2 assumes that the pressure distri- 
bution along the surface is sufficiently well known that an 
empirical correlation can be found to obtain 9p*/9x* | 
Conditions 3 and 4 state matching of the shear layer flow, 
expressed as an error function velocity profile 

u* = I U|[l - erf(C*)] (3.16) 


with the wake flow at the edge of the separated flow. 
Condition 5 states conservation of mass in the recirculation 
region, if the zero streamline is assumed to encompass the 
rear separation bubble. 

Further details and results of the incorporation of 
this wake flow model into the boundary layer approach can be 
found in Reference [31] and shall therefore not be reported 
here. The primary interest of this investigation is to 
study the flow through which landing and ascending aircraft 
would pass, i.e., the region somewhat above where the 
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boundary layer analysis described in the preceding is 
expected to provide meaningful results. 

IV. NUMERICAL SOLUTION 

The governing equations (3.3 through 3.5) were noh- 
dimensionalized by adopting some characteristic length L* 
and velocity U* from the flow field over the fence. The 
height/ h* , of the obstruction was chosen to be the 
characteristic length for this case. The characteristic 
velocity is somewhat arbitrary and is assigned a value equal 
to that of the undisturbed logarithmic velocity profile 
evaluated at three obstacle heights above the ground. The 
resulting equations in nondimensionalized form are then 
given by: 


3u . 9w _ 
3x Tz 


0 


-f 3u ^ — 3u 


■ 'lx i ^’1?) 


where 



X* 


z* 

X = 


Z 

■ 


u* 


w* 

u = 

U^' 

W 

U* 


oo 


00 

D = 

p* 


e:* 


PJU*2' 


and Re denotes the Reynolds number 


(3.17) 

(3.18) 


(3.19) 


(3.20) 


(3.21) 


I 



6a 


Re 


u*h* 

00 


(3.22) 


Expanding Equation 3.18, the equations to be solved are 
given by: 


Su . 3w 
dx 9z 


0 


(3.23) 


- 3u 


w 

d Z 3X 




3z: 


•d + e) 


3z 3zj 


(3.24) 


with the eddy viscosity model relating the turbulent motion 
to the mean flow variables: 


e 


(zo + z)^<^ 


du 

dz 


(3.25) 


These three equations together with the boundary and initial 
conditions form a closed set of nonlinear, parabolic, 
partial differential equations. They were approximated by 
an implicit finite difference scheme. After linearization 
of the inertia terms, the resulting tridiagonal matrix was 
solved by an elimination routine. A detailed description is 
given in Frost, et al. [21] . 

The lower boundary condition imposed a no-slip mean 
velocity at the wall. In addition surface roughness effects 
were incorporated into the numerical procedure by applying 
the logarithmic velocity distribution to the two points 
closest to the lower boundary, thereby implying a wall layer 
of height-in-variant shear stress and friction velocity u^^. 
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It was assumed that the value of u* is given by: 

(3.26) 


where u (3) is the velocity at the grid point z = 3Az. The 
velocity at grid point 2 can be calculated by the logarithmic 
law as: 

u (2) « (3.27)^ 

and u (1) will remain zero. 

The upper boundary condition specif ies -the velocity 
to be equal to the initial logarithmic profile at 10 fence 
heights above the obstruction. Thus the velocity profile, 
which merges with that of the logarithmic profile due to the ■ 
decaying pressure gradient, is matched with the undisturbed 
flow at the outer boundary. 

The velocity w in the vertical direction was obtained 
by integrating the continuity equation (3.23) using the 
above calculated velocity profile of u. 

The choice of mesh size for the problem considered 
was dictated primarily by the truncation error of the 
difference equations and, in some cases, by the stability 
criteria. A balance between the accuracy required in the 
solution and the computing time necessary to meet these con- ■ 
ditions had-^to be achieved. For these reasons plus some 
nimierical experimentati,on the following mesh sizes were 
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chosen. A constant grid size in the marching (x) direction 
of 0.1 the fence height h and a variable mesh size in the 
vertical (z) direction incorporated through a parabolic 
stretching function. Typical grid sizes ranged from 0.01 
near the wall to 0.35 of the obstacle height at the upper 
boundary. Computing times for a 20.0 x 10.0 flow region 
were approximately four minutes on an IBM-360-65 computer. 

V. RESULTS AND DISCUSSION 

Numerical solutions of the turbulent boundary layer 
equations have been carried out to assess the influence of a 
number of different parameters on the solution. The fol- 
lowing parametric effects are discussed in corresponding 
order: (1) the influence on the solution of the pressure 

variations taken along different nonseparating streamlines 

(f' _ / (2) the effect of the quadratic pressure decay function 

ns 

as compared to a more natural pressure decay based on the 
inviscid model of Chapter II, and (3) the change of the 
resulting flow field due to parametric variation of the 
approaching logarithmic wind profile. 

As to the first effect, calculations using the decay 
function given in Equation 3.9 were carried out for the 
t '? . f streamlines 

...f •; 

vsical location as calcula \ from 

the resulting pressure distribution 


--- 0.05; 0.1; 0 


lectively. The 

■ -vxS'.:'. 
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can be seen in Figure 2.11, page 42. Only small changes 

occur in the computed velocity profiles for the different 

solutions. It should be noted that the profiles presented 

here and in the following are as obtained in the assumed 

numerical coordinate system shown earlier in Figure 3.1b, 

page 52, i.e., the profiles of various x-stations are 

plotted relative to the same origin and are not 

shifted according to their actual location in the physical 

coordinate system. Comparing Figures 3.4 and 3.5 one finds 

that at the location of maximum velocity overshoot, at 

X * 6,0 there is only a 1.2% difference in velocity between 

the - 0.05 and the ip „ = 0.2 solution, with the velocity 

ns ns 

increasing toward the latter. It is apparent that the 
solutions are relatively insensitive to the choice of 
streamline which lends credence to the assumption that the 
far flow field is driven by the pressure gradient along the 
first streamline for which separation does not occur. 

The effect of the pressure decay based on the 
inviscid model (Equation 3.12) as compared to the quadratic 
pressure decay (Equations 3,8 and 3.9} was investigated 
eunong others for the - 0.05 streamline. The results are 
shown in Figures 3 . 5 and 3.6. It can be seen that the 
quadratic decay function dampens the overshoot and merges 
the velocity profile somewhat faster into the undisturbed 
logarithmic profile than the inviscid decay function. The 
slower decay of the inviscid model seems more realistic as 
it distributes the disturbance better in the vertical 




Heigh 



0‘0 0>a 0>«i O'B 0*3 i*0 i-2 1-4 1-S 

Velocity u 


Figure 3.5. Velocity profiles for = 0.05 and quadratic 
pressure decay function. 
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direction. However, the difference between the two is small 
and only a comparison with experimental data can show which 
is the better decay function. 

The combinations of parameters investigated in the 
third group to ascertain the influence of the approaching 
wind on the solution are tabulated in Table 3.1 and graph- 
ically presented in Figure 3.7, which is similar to the 
previously discussed Figure 2.5, page 33. The respective 
data points (Cases 1-4) were selected in the following way: 
four different surface roughness parameters zo ranging from 
Zo ~ 0.005 to Zo = 0.05 were chosen. With Equation 3.26 and 
the earlier assumption that the reference velocity is that 
at three obstacle heights above level terrain one can deter- 
mine the corresponding dimensionless friction velocities u* 
from: 


u; 

U* 

oo 


in 


3.0 + z 0 1 
Zo J 


(3.28) 


The resulting input parameters Cp^ introduced into the 
inviscid model, determining the respective pressure vari- 
ation, is then found from Equation 2,65. Table 3.2 gives a 
complete list of the various empirical input parameters for 
the inviscid model and the corresponding calculation 


parameters. It should be noted, that y^, C 


■max 


and x„ were 


kept constant in all four cases. The final results are 
given in Figures 3.8 through 3,11. Comparing the various 
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TABLE 3.1 

PARAMETERS FOR THE DIFFERENT DATA POINTS INVESTIGATED 
Case Z 0 u* Cpj^ 


Boundary Layer Approach 


1 

0.005 

0.0625 

-0.7920 

2 

0.010 

0.0700 

-0.7693 

3 

0.020 

0.0792 

-0.6978 

4 

0.050 

0.0973 

-0.4805 


Inviscid Test Case 


Set 2 


0.005 


0.0647 


-0.8500 
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TABLE 3.2 

EMPIRICAL INPUT AND RESULTING CALCULATION PARAMETERS FOR 

RESPECTIVE DATA POINTS 


Case 

Cpb 

Uo 

cosh $ 

cosh (S 


1 

-0.7920 

0.70711 

1.3348 

10.705 

1.00830 

2 

-0.7693 

0.70711 

1.3436 

10.697 

1.00170 

3 

-0.6978 

0.70711 

1.3733 

10.669 

0.98091 

4 

-0.4805 

0.70711 

1.4827 

10.566 

0.91163 
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Figure 3.8. Velocity profiles for case 1. 
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Figure 3.10. Velocity profiles for case 3. 
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velocity profiles for the two extreme cases 1 and 4, one 
finds that for increasing surface roughness Zg and friction 
velocity u^^ there is a decrease of velocity overshoot of 
about 7% at almost all x-stations over the recirculation 
region downstream of the obstruction. Whereas the roughness 
effect, i.e., decreasing velocity for increasing roughness 
is dominant in the inner flow region, the influence of the 
friction velocity, i.e., increasing velocity for increasing 
u^ dominates in the outer flow field. This tendency appears 
realistic as the fuller approach velocity profile with 
smaller Zo, which carries more mass and higher momentum near 
the wall can be expected to have the larger overshoot. 

Plotting the profile of the velocity increase due to 
the presence of the obstacle near the highest point of the 
recirculating region, i.e., at x = 6o0 shown in Figure 3.12 
one observes the same trend: Case 1, with the fattest 

initial profile due to its lowest surface roughness, shows 
the largest velocity increase. 

This suggests that a high rise building, for example, 
located in a downtown area of large buildings, i.e., high 
surface roughness, will not experience as significant an 
overshoot as the same building in a resort or residential 
area surrounded by natural terrain or small houses having 
low surface roughness. 

With regard to the overshoot in the velocity pro- 
files, one would generally expect it to be somewhat smaller 
in magnitude and to be distributed over a wider vertical 
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range r as the full scale experimental data of Frost, et al. 
[32] and partially the wind tunnel data of Good and 
Joubert [27] indicate. A comparison with these experiments 
is given in Figures 3.13 and 3.14. It should be pointed out 
that the former tests were conducted for a rectangular block 
building, whereas, the latter were for a fence. In order to 
make the respective velocity distributions independent of 
the specific parameters of the approach velocity and to 
facilitate a comparison, all profiles are referenced to 
their undisturbed upstream profile. In Figure 3.13 the 
boundary layer profiles of the previous Figure 3.5 (qua- 
dratic decay), page 65, and Figure 3.6 (inviscid decay), 
page 66, are plotted together with a typical profile obtained 
by Reference [32] in an atmospheric boundary layer flowing 
over a block building. The horizontal location of these 
profiles is one obstacle height downstream of the obstruction 
face. The stronger and more confined overshoot region pre- 
dicted by the boundary layer model probably results from the 
fact that the wake is treated as a solid body. Because it 
is intended to match the present boundary layer model with 
the proposed integral technique for the recirculation region 
later, this was believed to be a reasonable temporary 
assiamption. It implied, however, that not only the velocity 
decays to zero but also that the mixing length reduces to 
£ * K • zo along the separating streamline, leading to a 
relatively small eddy viscosity and low shear in this region, 
not at all representing the highly turbulent free shear 
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Figure 3.14. Comparison of velocity profiles at x = 2 
obtained from boundary layer model and wind tunnel 
experiments [27] . 
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layer. The result is the larger virtually laminar overshoot 
in the vicinity of the separating streamline. High shear 
would not only reduce the overshoot in magnitude but also 
spread it out further in the vertical direction. The mass . 
and momentum diffusion from the shear layer along the 
separating streamline into the recirculation region cannot , 
be taken adequately into account through the solid surface 
assumption and thus lead to the increased overshoot. 

On the other hand, comparison with the small scale 
tests of X27] in Figure 3.14 shows a measured overshoot 
which is even larger than that predicted by the theoretical 
model, however, less rapidly decaying in the vertical 
direction. 

It should be noted here that the two experiments 
differ in an important parameter. While the first tests 
(32] were conducted in an atmospheric boundary layer where 
the ratio of obstacle height h* to boundary layer thickness 
6 * is very small, i-e.. 


h* 

T*- 


<<• 


1 


the latter tests were conducted in a smooth flat plate type 
boundary layer with 



2.5 


and hence should, because of the increased mass displacement 
near the obstruction, produce a larger overshoqt. Taking 
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this into consideration, the predicted velocity profiles for 
atmospheric flows should tend towards the first experiments. 
If, in addition, it is realized that these tests are not 
strictly two-dimensional because of the finite width of the 
block building, but have three-dimensional character 
allowing flow around the sides, thereby diminishing the 
overshoot over the top of the building, the comparison with 
the calculated velocity profiles may look more reasonable. 

Concluding this chapter, one might say that because 
of the solid body wake assumption and the related simplif i- • 
cations, e.g., in the mixing length prescription, the 
boundary layer analysis yields results which show overshoot 
regions in the velocity profiles that are somewhat higher 
than expected. As there are not sufficient experimental 
data available for direct comparison, only qualitative 
conclusions can be drawn, following the trends pointed out 
in the foregoing discussion. 


CHAPTER IV 

ANALYSIS OF ATMOSPHERIC FLOW OVER A BLUFF 
SURFACE OBSTRUCTION BY THE TURBULENT 
NAVIER-STOKES EQUATIONS 

The complete two-dimensional equations of motion are 
applied to an atmospheric flow over a forward facing stepr 
as shown in Figure 4.1. Analogous to the previous approach 
the atmospheric motion far upstream of the obstruction is 
described by the logarithmic profile given in Equations 3.1 
and 3.2. Again the atmosphere is assumed to be neutrally 
stable and Coriolis effects are assumed negligible due to 
the small scale of motion under consideration. 


I . GOVERNING EQUATIONS 


With the above assumptions the governing turbulent 
mean-flow equations for steady incompressible flow can be 
written as follows. 

Momentum equation in x*-direction 


u* 



+ w* 


au* _ _ 1 ap* , a 

r * 

^eff 


2 

az* p* ax* ax* 

p* 

Mo 

1 3X*J 


+ 


r * ^ ' 

a ^eff f au* aw*1 

az* p* [az* ax*J 


(4.1) 
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Figure 4.1. 


Description of flow region considered. 
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momentum equation in z* -direct ion 


u* 



+ w* 



1 3P* 9 

p* 3z* 3z* 

3 f^efffSu* ^ 3w*ll 



and the continuity equation 


3u* 3w* _ 

3x* 3z* ^ 


(4.2) 


(4.3) 


where overbars again denote the ensemble average. Analogous 
to the turbulent boundary layer concept, turbulence is 
incorporated through an "effective viscosity" 

y*ff = y* + y* (4.4) 

which is composed of the molecular viscosity y* and a turbu- 
lent viscosity y*. Unlike y*, the turbulent viscosity y* is 
not a property of the fluid. Its value will vary from point 
to point in the flow, being largely determined by the 
structure of the turbulence at the respective location. 
Several methods, discussed in detail in the next paragraph, 
are used to express the turbulent viscosity in terms of 
known or calculable flow quantities. To eliminate the 
pressure from the governing equations (4.1 and 4.2) the 
stream function ip* is introduced 

u* = Ml 
“ p* 3z* 


(4.5) 
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w* = 


1 3if>* 

F*” 3x* 


together with the vorticity w* 


(4.6) 


3w* 3u* 

3x^ 


3z^ 


(4.7) 


Differentiating the x*-momentuin equation with respect to z* , 
the z*-momentum equation with respect to x* and subtracting 
one from the other, one obtains after some rearrangements 



(4.8) 
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3x*5| 

^eff TP^I 





< * 


32 


'n* -Kl 
J*^eff ^x^J 


j 


(4.9) 


Equation 4.8 is generally known as the vorticity transport 
equation. It reduces the two momentum equations (4.1 and 
4.2) in the three variables u* , w* , p* , to one equation in 
the two variables w* and ip*. 

A second equation can be derived from Equation 4.7 
with Equation 4 . 6 


V^\p* = -p*w* 


(4.10) 
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known as the stream function equation. Gosman, et al. [20] 
have shown that these new governing equations (4.8 and 4.10) 
can be expressed in the common form of an elliptical partial 
differential equation, suitable for simultaneous numerical 
integration , i . e . , 



- 5|r(b + a - 0 (4-11) 


Here (j) is the respective dependent variable ifi* or co* and a, 
bf c and d are functions depending on the variable under 
consideration . 

After rearranging some terms. Equation 4.11 can be 
rewritten, using tensor notation, as 


ap*u^ 


9 (|) _ 9 

9x| 9x^ 


b 




d 


(4.12) 


Table 4.1 lists what the functions "a" through "d" must be. 


A. 

( 0 * 

Ip* 


TABLE 4.1 

COEFFICIENT FUNCTIONS OF EQUATION 4.11 OR 4.12 

a b c d 


1 


‘eff 


-S* 

(d 

“( 0 * 


0 


1 

1/p* 


1 
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Gosman/ et al. [20] developed an algorithm for solving 
Equation 4.11 which is used in this study. 

II. TURBULENCE MODELS 


Following the turbulent viscosity concept, there 
remains the task of formulating y*, that is relating it to 
some known quantities of the mean flow. These resulting 
auxiliary relationships can either be simple algebraic 
expressions or more sophisticated differential equations. 
Together with the governing equations (4.8 and 4.10) and 
appropriate boundary conditions they form a closed set of 
equations describing the flow for the problem of interest. 


Mixing Length Model 

Among the models which employ algebraic relations 
for y* is Prandtl's mixing-length hypothesis (PML) , already 
known from the boundary layer approach. The hypothesis is 
that the turbulent viscosity is equal to the local product 
of the density, the magnitude of the mean rate of strain and 
the square of a characteristic length scale of the turbulent 
motion, i.e., the mixing length I*. 


y 


* = 
t 






3u* 

3z* 


3w* 


3x* 


1/2 


(4.13) 


The mixing length must be prescribed algebraically for the 
whole flow field. In the present investigation its distri- 
bution is made dependent on the shortest distance from the 
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wall as indicated in Pigure 4.2. 


l*(x*,z*) = 


ic(z* + z*) 


K ( |x* I + z*) 


ic(z* - h* + z*) 


for X* < 0; 
z* < |x*| 

for X* < 0; (4.14) 

I X* I < z* < |x* I + h* 

for all X*; 

Z* > |x*l 


Turbulence Kinetic Energy Model 

We now turn to a model of turbulence where the 
determination of y* requires the solution of a differential 
equation for one property of turbulence. The model was 
first suggested by Kolmogorov (33] and Prandtl [34] and 
differs from the previous mixing length model by the assump- 
tion that y* is dependent on the level of turbulence of the 
fluid and a length scale. The level of turbulence is 
characterized by the mean kinetic energy of the velocity 
fluctuations, defined as 


k* = j(u* * 2 + V* * 2 + w* ' ^) 


(4ol5) 


where u* ' , v* ' and w* ' are the fluctuating parts of the 
velocity components. The quantity k* is called the turbu- 
lence kinetic energy (TKE) . It is related to the effective 
viscosity by 

y|ff = P* • • I* • 


(4.16) 
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where is an empirical constant or function to be dis- 
cussed later, k* is determined from a differential equation 
of the same form as Equation 4>12, containing convection, 
diffusion and source terms. It may be shown after some 
lengthy algebra that such an equation can be derived from 
Equations 4.1 through 4.3. For a detailed derivation see 
Wolfshtein [35] . This transport equation for the turbulence 
kinetic energy is 


p*u^ 


3k* 

3x* 

3 


3x^[ k,eff 3x^J 


^ ^k 


(4.17) 


where r* is an exchange coefficient for the turbulence 

kinetic energy defined in terms of a Schmidt number based on 
the effective viscosity 


f* 

k,ef f 


K±f 

’k,eff 


(4.18) 


The last term in Equation 4,17 represents "sources” and 
"sinks" of turbulence; it consists basically of two parts, 
one accounting for the rate of generation of turbulent 
kinetic energy by the turbulent shear stress and another 
representing the energy dissipation by viscous action, 
deducible from dimensional analysis 









k*3/2 

Z* 


(4.19) 


Cp is another function to be determined empirically 
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In this model the length scale 1 * is still pre- 
scribed algebraically and is assumed to be equal to the 
mixing length distribution of Equation 4.14. 

The determination of the functions and still 
poses some problems especially in cases where there are no 
experimental results available. In terms of a turbulent 
Reynolds number 


R 


t 


p* » « A* 

y* 


(4.20) 


one can express the functions as 



i.e., when is large we can expect the molecular viscosity 
to have negligible effect on the transport process and the 
functions take on constant values. However, when is very 
small and turbulence effects are negligible, the functions 
tend to 

«t 

The "universality" of the constants Cy^ and Cjj^, for 
which many users of the model had hoped, could not be 
achieved; the values were found to be quite different 
depending on the experiment under consideration. Wolfshtein 
135] , modeling a jet impinging normally on a wall, uses 
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Cy, = 0.22; Cdj, = 0.416 

Ng and Spalding [36] use 

Cy ^ = 1.0; Cp ^ = 0.1 

for their calculations of boundary layers near walls, as do 
Rodi and Spalding [37] in their modeling of free turbulent 
flows. Launder, et al. [44], however, employ 

Cy u = 0.09; Cp = 1.0 

in their niamerical solution for free turbulent shear flows. 

In view of the fact that there is only insufficient 
experimental data available for the present investigation, 
no attempt shall be made to propose new values for the two 
constants in question. Instead, the influence of their 
variation on the flow field shall be studied. 

Two^Equation Model of Turbulence 

The next logical step from the turbulence kinetic 
energy model described in the previous paragraph is to 
remove the uncertainty from the length scale distribution 
especially in recirculating flow regions and to calculate i* 
rather than to prescribe it algebraically. Rotta [38] was 
the first to propose a differential equation for the length 
scale, deriving it from the Navier-Stokes equations. 

Gosman, et al. [20] reformulated the equation into the 
common form of their transport equation, thus, together with 
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the k* equation they arrived at a two-equation model of 
turbulence o The equation for is 


P*u5 


r + a* 

3x*l^il,eff 3x*J 
3 '■ 3 ^ 


(4.22) 


with 



3x? 3x? k» 
^.1 ^ 


Cg + p*k*^/^ 


(4.23) 


being the ” source” term of the turbulent length scale. It 
consists of a positive contribution representing the rate Of 
growth of Z* as a result of the dissipation of -energy, 
especially from smaller eddies and a negative contribution 
accounting for the tendency of the shear stress to reduce Z* 
by rupturing the large eddies. Cg and Cg are again ftinc- 
tions of having similar behavior as and Cg in Equation 
4.21, x.e.. 



(4.24) 


C 3 (j and are constants to be determined empirically. 
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III o NUMERICAL SOLUTION 
Numerical Procedure 

The mathematical problem posed in the previous 
chapter has been solved by utilizing and extending the 
numerical procedure of Gosman., et al. [20] . For a detailed 
description the reader is referred to this reference. Only 
a brief outline is given here. 

The governing differential equations (4.8, 4.10, 

4.17, 4.22 or 4.12 in general) are replaced by algebraic 
finite difference equations which are obtained by integration 
over finite areas rather than Taylor series expansion, 
assuring a broader range of applicability especially in non- 
rectangular coordinate systems. The integration of the 
convection terms employs "upwind differencing," a one-sided, 
rather than centered space differencing, where the scheme is 
backward when the velocity is positive and forward when it 
is negative. This formulation of the first order terms 
gives greater numerical stability than can be obtained with 
central differences. The remaining diffusion and source 
terms, however, are expressed in a weighted central differ- 
ence form. 

Because of the nonlinear character of the resulting 
finite difference equations they are solved by an iterative, 
successive substitution technique. 
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Physical and Numerical Coordinate System 

The physical coordinate system was chosen such that 
the origin was located at the lower step corner, with the 
positive x*-axis pointing in downstream direction parallel 
to the wall and the z*-axis directed normally to it aligned 
with the front face of the step.*. (See Figure 4.1, page 83.) 
In this coordinate system the flow regime considered extends 
ten step heights in the upstream and downstream directions 
and nine step heights in the vertical direction. The origin 
of the numerical coordinates was situated at the lower left 
corner of the flow field (Figure 4.3) with I indexing in the 
x*-direction and. J indexing in the 2 *-direction. In the 
iteration process the field was swept from left to right 
beginning at the wall and proceeding in the increasing J- 
direction. 

The distribution of the grid points is shown in 
Figure 4.4. As indicated, a variable mesh was used which 
gradually decreased in size near the wall and in the 
vicinity of the step* 

Boundary Conditions 

Depending on the turbulence model under conside- 
ration the number of differential equations to be solved 
ranges from two equations in oi* and ip* to four equations in 
CO*, i|)*, k* and A*. The number of boundary conditions 
required for the respective case then changes accordingly. 
Due to the elliptical nature of the flow problem, these 
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boundary conditions are prescribed along the entire boundary 
of the flow regime/ i«e./ following Figure 4.1/ page 83/ 
along the inlet/ the outlet/ the upper and lower boundaries 
(wall) . All conditions are either of Dirichlet or Neumann 
type. 


Inflow . At the inflow a logarithmic velocity pro- 
file of the form 


u* = 



z* + z* 

Z * 


(4.25) 


is assumed. The boundary condition can then be deter- 
mined by integrating the velocity profile over the inlet 
height 


Ip* (z*) 



z* 


/ 


fn 


z* + 


z 


7 

0 


* 

Zo 


dz* 


(4.26) 


which yields 


uJ 


ip* (z*) = p* 


(z* + z*) Jin 


•X. » * 

Z* + Zq 


(z* + z*) 


J 


+ C (4.27) 


where C can be determined such that = 0 at the wall. 
Then 


u2 


ij^*(z*) = p’ 


(z* + z*) 


Jin 


Z* + Z( 


+ Z* 


(4.28) 


The condition for the vorticity to* is 
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( 0 * (z*) 


_ 9w* 3u* 

“ Sx* ” 9z^ 


(4.29) 


where the last term can be calculated from the ip*- 
distrlbution as 


9u* ^ 1 ^ ^ 1 

3z* ^ 9z'*"2' K z* + z* 

Cf 


(4.30) 


The remaining term of the inflow vorticity was 
allowed to develop as part of the solution by approximating 


3 ^w* _ 
3x*2 


0 


(4.31) 


and setting 


3w* 


3x* 


1,J 


3w* 
3x* I 


2,J 



3 

3x*2 


2,J 


(4.32) 


The employed boundary condition for the turbulent kinetic 
energy k* was derived from the constant shear assumption 
underlying the derivation of the logarithmic law of the 
wall. 

T*(z*) = p*u*^ = const. (4.33) 


In terms of the Prandtl-Koimogorov formulation (Equation 
4.16) one can write 


T* (z*) 


p*k*^/^ic*C 


M 0 


3u* 
3 z* 


(4.34) 


Equating the two and substituting from Equation 4.14 and 
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the velocity gradient from Equation 4.30 one obtains the 
boundary conditions for k* at the inlet as 


k* = 



(4.35) 


The corresponding condition on the turbulence length scale 
is then 


I * - k(z* + z*) (4.36) 

Outflow . Two sets of outlet boundary conditions 
were used. The first one assxnnes that the outlet location 
is sufficiently far downstream of the step that an undis- 
turbed logarithmic velocity profile has developed again. 


U. . in * I* - h* 


(4.37) 


can then again be determined by integrating u* 


u 


ij;*(z*) = 


*out 


(z* + z* 


- h*) 


Jin 


:* + z* 


- h* 




+ z* 


(4.38) 


where is calculated from the conservation of mass 
through inlet and outlet/ as no mass will be assumed to 
cross the other boundaries. 


V* = 
in out 


(4.39) 


with 
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H* 


VI' 

xn 


-I 


u*dz* 


Vt 

xn 




(H* + z*) An 


H* 


* 

+ Z Q 
w * 


- H* 


( 4 . 40 ) 


and 




z* 

* 0 


(H* - h*)j 

( 4 . 41 ) 


the friction velocity at the outlet is 


u 




*out 


(H* + z* - h*) An 


H* + z!s - h* 


)- 


( 4 . 42 ) 


(H* - h*) 


Analogue to the inlet, the first term of the vorticity 
boundary condition is determined from the ip* distribution by 


3u* _ 1 
3z* ^ 9 z *2 


( 4 . 43 ) 


while the other co*-term is calculated with the afore- 
mentioned assumption 


3^w* 


0 


( 4 . 44 ) 


1 • 3 o f 
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p* 3x^2 


( 4 . 45 ) 
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The remaining conditions on k* and I* were similar to 
Equations 4.35 emd 4.36 


k* 


^*out l 

Cyo J 


2 


(4.46) 


and 


I* • K(z* + z* - h*) (4.47) 

The second set of outlet boundary conditions consisted of 
less restrictive formulations for w* and ip*. 

The assvunption for ip* was that 


d^w* _ 
3x*2 * 


0 


(4.48) 


or 


a*** 

d3^ 


IN,J 


B^ib* 

dx*^ 


IN-1, J 


’^'IN - ’^IN-2 



(Ax*)^ 


which yields 


’^IN,J 


* 1'™- 


IN-2 


(4.49) 


For the vorticity it is assumed that its gradient in the 
flow direction vanishes 




0 


(4.50) 


The turbulent kinetic energy k was allowed to decay to the 
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initial free stream turbulence level 


k* = 



2 


(4.51) 


Upper boundary . The location of the upper boundary 
was assumed far enough out that velocity deflections caused 
by the step were negligibly small. Consequently the stream- 
line condition is 


Ip* as const. (4.52) 

i.e., no flow is crossing the upper boundary. 

The vorticity condition imposed was that of a 
vanishing gradient 


3z* 


0 


(4.53) 


which was also required for the turbulence kinetic energy , 
however, in the horizontal or x*-direction 

® 0 or k* = const. (4.54) 

The length scale was prescribed as given by Equation 4.14. 

Wall boundary . As the flow is parallel to the wall 
the condition for the stream function must be 

Ip* » const. (4.55) 

along the entire wall boundary. 
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The vorticity boundary condition is more problem- 
atic, as it essentially drives the flow. It can be derived 
from a Taylor-series expansion of the stream function around 
a near wall point (NP) , An away from the wall, in terms of 
the wall point (P) conditions 



+ 



An* + 
P 


1 

2 3n*2 


(An*)^ 

P 


3,1,* 


+ i i!it_ 

6 3n*3 


(An*) 


+ H.O.T. 


(4.56) 


By the no-slip condition: 


1^=0 


(4.57) 


3^ij^ * _ 
3n*^ 


p*w* 


(4.58) 


Combining the last two expressions, one gets 


3^ii<* ^ __3_ 
3n*3 3n* 


- 

[3n*2j 


* 3o)* 


(4.59) 


Substituting these into Equation 4.56 and solving for w* 


wi = ~ 
p 


- »g) 

(An*) • p* 


3u*| An* 




(4.60) 


The vorticity gradient at the wall was approximated by 



(I). 


NP 


- (O’ 


An* 


(4.61) 


yielding the second order vorticity formulation 
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OJ? 


(An*)^ 



(4.62) 


The turbulence kinetic energy specification along 
the wall poses some problems, k* = 0 is incompatible with 
the inlet and outlet conditions of 





and 




u 


*out 


'Po 


(4.63) 


(4.64) 


It was therefore assumed that k* obeyed Equation 4.63 also 
along the wall, however, with uj varying as 



wall 


i* 


3u* 

dz* 


wall 


(4.65) 


implying a logarithmic velocity profile from the wall to the 
first interior node. Furthermore, at the wall, 

“Jlwail “ K‘z*»(-u)*) (4.66) 


or 


^wall 


_ Tic * z**g)* 

i Cy 0 


(4.67) 


The length scale was prescribed by 
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I* = K«z* (4.68) 

Special treatment at step corner . The stream 
function at the upper corner presents no problem. Like the 
rest of the wall = 0. Similarly, k* and Z* are uniquely 
prescribed by the lower boundary condition. But there are 
several alternatives for the evaluation of w*. Referring to 
Figure 4.5 one can apply Equation 4.62 either to the up- 
stream side (face) of the step or to the downstream wall 
obtaining respectively: 


0 ) 


* SB 


= 10 


B 




- K CO, 


p* (Ax*) 


2 '"W 


(4.69) 


and 



^I’n . 1 

p*(iz*)2 2 


( 0 , 


N 


(4.70) 


There are also other possibilities. Seven different methods 
are given in Roache [39] , and most of these are investigated 
in the present study; they are listed in Figure 4.5. Method 
(1) represents an attempt to force separation at the corner, 
assuming that the vorticity vanishes at a separation point. 
Method (2) is based on the idea that, since separation 
occurs tangentially to the upstream wall, upstream wall 
evaluation should be used. Method (3) is an attempt to 
average out both values, while method (4) is derived from 
adding both values in a first order formulation. Method (5) 




( 1 ) 0 )* = 0 


(4) w* = “2 (ip* + t/;*)/An*^ • p* 


(2) 0)* = CO* 

(3) ~ + w^)/2 



Figure 4.5. Investigated boundary conditions for upper step 
corner point. 


o 


as * 








arises from the argument that no continuity of to* can be 
expected for the geometric singularity of the corner. 

When applied to the laminar flow case^ obtained by 
setting the turbulent viscosity equal to zero, all methods 
functioned well and enforced separation from the corner. 

The second formulation, however, was the most effective in 
producing a well developed realistic separation region. For 
this reason it was used in all subsequent calculations. 

Accuracy, Convergence and Economy 

Carrying out the numerical solution to the above set 
of differential equations, a balance is required between the 
convergence and accuracy and the amount of computing time 
necessary to meet these conditions. One of the most 
important factors affecting this balance is the grid spacing. 
Therefore, the following sections will discuss the effects 
of grid dimensions on accuracy and convergence. 

Influence of grid size on accuracy . Accuracy is the 
deviation of the numerical solution from the exact solution. 
Unfortunately, this exact solution is unknown in most cases, 
thus the accuracy is not easily determined. One possibility 
to overcome this difficulty is to test the numerical pro- 
cedure with a simple problem for which the exact solution is 
known. This allows an estimate of the general accuracy of 
the procedure. The optimum grid size, however, cannot be 
transferred to another problem, as it strongly depends on 
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the physical case under consideration, i.eo, the size and 
location of the occurring gradients » 

Another way is to investigate the dependence of the 
respective solution on increment size and find the grid size 
for which a further decrease brings no further or at most 
minor improvements in the solution. The grid size for the 
present analysis was essentially determined by this approach 
with some additional compromise toward economy. The fol- 
lowing discussion shall compare some results obtained from 
the grid size actually used for all subsequent calculations 
with results obtained from a grid with twice as fine a mesh. 
Such a grid increases the computing time by a factor of 
four. All calculations were carried out for a surface 
roughness z^ = 0»45 [m] and a friction velocity of u{ = 

0.75 [m/sec]. Figures 4.6 through 4,11 show the vertical 
vorticity and stream function distributions at three x*- 
stations for the two different grid systems. In the first 
four figures for the x* = 6.0 h* and x* = 0=0 locations 
there are almost no differences in the respective oi*- 
and Tj;*-distributions and consequently no improvement through 
the finer mesh. The situation is slightly different a short 
distance behind the step corner at x* = 0.8 h*. The stream 
function, the better behaved function of the two, is again 
almost identical in both cases. The vorticity-distribution, 
however, shows differences. Even the finer mesh does not 
quite predict the wall vorticity correctly, although it 
shows the correct trend, indicated by the dotted line in 
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Figure 4,6. Stream function distribution at x 
different grid sizes. 
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Figure 4.7. Vorticity distribution at x = -6.0 for 
different grid sizes. 
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Figure 4.8. Stream function distribution at x = 0.0 for 
different grid sizes. 
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Figure 4,10. Stream function distribution at x = 0.8 for 
different grid sizes. 
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Figure 4 <,11. Because of the expected flow separation at the 
top step corner and the resulting recirculation region 
behind it, the vorticity should reverse its sign in the 
vicinity of the wall as one approaches the wall in the 
vertical direction from the inside of the flow field. This 
can best be demonstrated by looking at some results from a 
laminar flow case obtained by setting the turbulent vis- 
cosity to zero. Here the separation region is bigger and 
thus more grid points fall into this region. Looking at 
stream function distributions for two different x*-stations 
downstream of the step corner in Figures 4,12 and 4.14, and 
a vorticity distribution in Figure 4,13 corresponding to the 
first x*-station, it is seen that at x* = 1.6 h* the coarser 
grid does not give negative streamline values near the wall, 
which one would expect because of the separation from the 
corner. The finer grid, on the other hand, produces the 
expected negative i|^*-values. At x* = 3,0 h* the separation 
region is large enough that even with the coarse grid 
negative i/;*-values can be obtained near the wall. Conse- 
quently the stream ‘function plots for the two grids look as 
shown in Figure 4,15, where (a) is the plot for the coarse 
and (b) the plot for the fine grid. The two vorticity 
distributions for the x* = 1.6 h* stations are shown in 
Figure 4,13. They demonstrate how inside the separation 
region the vorticity changes sign in the vicinity of the 
wall even for the coarse grid. Because of the stronger 
developed separation region in the laminar case, i.e.. 
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Figure 4.12. Near wall stream function distribution for 
laminar flow at x « 1.6. 
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Figure 4.13. Near wall vorticity distribution for laminar 
flow at X = 1.6. 
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Figure 4.14. Near wall stream function distribution for 
laminar flow at x » 3.0. 
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greater physical extent, vorticity wall predictions are 
better than in the turbulent case where the separation 
region is considerably smaller. 

The above fine grid calculations were obtained from 
solving a small flow regime behind the step with fixed 
boundary conditions, interpolated from coarse grid calcu- 
lations, except for the vorticity boundary condition at the 
wall, which was allowed to develop as a new solution. 

It should be pointed out that the oj*-boundary con- 
dition (Equation 4.62) is based on the linear vorticity 
distribution. However, as can be seen from previous figures, 
the vorticity is linear only in the very close proximity of 
the wall. When a fine enough mesh is used, the point near 
the wall is in this linear vorticity region and the wall 
vorticity predictions are good. In either case, the vor- 
ticity away from the wall is not very sensitive to mesh size 
and is hardly influenced by the local inaccuracies at the 
wall inside the turbulent separation region. 

To summairize the findings of the foregoing investi- 
gation, one might say that except for a small region in the 
immediate vicinity of the wall at a short downstream section 
behind the step corner the coarse grid produces reasonable 
results, which can only be modestly improved by the fine 
grid calculations requiring four times as much computing 
time. The calculations show further that inaccuracies in 
the vorticity wall predictions inside the turbulent sepa- 
ration region on top of the step, do not propagate far 
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Inside the flow field but remain limited to a small area 
near the wall. The external flow field remains largely 
unaffected. 

If one wishes to improve the predictions in the near 
wall separation region, however, one can choose from several 
possibilities. The first has already been mentioned during 
the discussion, i.e., decrease of the grid size near the 
wall until the variation of all computed quantities changes 
approximately linearly between adjacent mesh points. 
Unfortunately, this is not always easily accomplished. 
Besides sizeable increases in computing time, convergence 
problems can easily occur when unsuitable grid distributions 
are used. More information on this subject is contained in 
the next section dealing with convergence. 

A second possibility is to leave the grid unaltered, 
but abandon the assumption of linear variation of properties 
near the wall, putting in its place some information about 
the way in which the properties actually vary in the inter- 
val in question. The relations containing this information 
are commonly referred to as "wall functions . " Suitable 
formulae for turbulent flows near smooth walls with zero 
pressure gradient were derived from Couette flows by 
Patankar and Spalding [40, 41] and Wolfshtein [35] . While 
the former base their functions on the Prandtl mixing-length 
concept, the latter makes use of the Prandtl-Kolmogorov 
hypothesis. In both cases, some of the important constants 
and functions needed for the completion of the system of 
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equations were deduced by reference to experimental data. 

For rough walls ^ no wall functions have been derived 
yet. A paper by P. A. Taylor and Y. Delage [42] is a first 
approach in this direction. For their computation of 
atmospheric boundary layers with zero pressure gradient over 
rough terrain they assume a constant flux wall layer having 
a logarithmic velocity profile for the calculation of the 
first interior grid point. The formulation of the wall 
boundary condition for the turbulence kinetic energy 
(Equation 4.67} is based on this assumption. 

The provision of a comprehensive set of wall func- 
tions valid for most situations of practical interest is one 
of the prime tasks of current research in computational 
fluid dynamics dealing with turbulence. 

A third alternative » much less sophisticated, 
follows the basic wall function concept to account for the 
nonlinear behavior near the wail. Taking advantage of the 
fact that in the present problem w* and \p* are fairly 
independent of grid size inside the flow field it determines 
the location of the separating streamline, i.e., the second 
zero of Ip*, by extrapolation from the interior of the flow 
field. Figure 4.16(b) shows a typical turbulent flow ip*- 
distribution through the separation region as obtained with 
the coarser grid. With the assumption of linear variation 
the distribution between grid points follows the solid line. 
This is a good approximation except for the small region 
near the wall, where the estimated correct distribution is 
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represented by the dotted line. The second zero of ip*, the 
separation streamline # will always be neglected by the 
linear approximation if it falls between the first two grid- 
points. Because of the almost linear^ slightly parabolic 
behavior of the stream function closer to the wall, as seen 
in Figures 4,6, 4,8, and 4. '10, pages 110, 112, and 114, a 
parabolic extrapolation was used for the location of the 
second zero of the stream function, starting from the first 
two interior ip values towards the wall. A typical result is 
presented in Figure 4,17 which shows streamline plots for 
the same flow case with (4.17 (b)J and without (4.17(a)) 
extrapolation of the separating streamline. 

Factors affecting convergence . While accuracy 
greatly depends on mesh size, the convergence of the 
iterative solution procedure heavily depends on the mesh 
size variation. Recommendations from an accuracy viewpoint, 
e.g., variable grid size near walls, have to be applied with 
caution when laying out a suitable mesh. It has been 
experienced by Gosman, et al. (20] and also in this study 
that near walls nonuniform grid spacing between grid lines 
parallel to the wall may cause divergence due to the 
coupling of the vorticity and the stream function equation 
through the vorticity boundary condition (Equation 4.62). 

The suggested [20] remedies listed below have been found 
quite effective, 

1. Near the ^&all the ratio of consecutive intervals 
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between the nodes normal to the wall should be 
kept as close to unity (i.e., uniform) as 
possible or otherwise below 1.5. 

2o The vorticity boundary condition (Equation 4.62) 
should not be used explicitly at the wall, but 
be incorporated into the general substitution 
formula for the iterative solution (implicit 
formulation for vorticity) . 

In the present investigation this last condition was not 
necessary to assure convergence. 

Another source of divergence is that even inside the 
flow field large variations of the coefficients in the sub- 
stitution formula may occur. This is true especially in the 
turbulent kinetic energy solution cycle for the source tern 
Sj^ (Equation 4,19) c As shown by [20], the substitution 
formula can in this case be rearranged through simple alge- 
braic manipulations such that variations in a modified 
source term stay small. It was found in the present study 
that this approach may also be used for the substitution 
formula of the turbulent length scale if its source term 
(Equation 4.23) should cause divergence, 

A more commonly employed remedy against divergence 
of the iteration process, although more time constiming, is 
known as under-relaxation. Compared to the method just 
discuss id, it is more easily applied in the various cases, 
but has the disadvantage of slowing down the iteration pro- 
cess, thus increasing the computing time, depending on the 
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degree of under-relaxation , For the present problem there 
was no need for under -relaxation. On the contrary, the 
stream function, related to the vorticity by a Poisson type 
equation (Equation 4.10), could be over-relaxed to speed* up 
the solution. 

Termination of computation . The computation was 
assumed to have converged to a sufficiently exact solution 
if the difference of the dependent variable (j) at point P 
between successive iterations became small, i.e.. 


.N .N-1 

<|)p - <|)p 




N 


< e 


max. 


0.01 - 0.05 


(4.71) 


The superscript N denotes the nth iteration. (})p rather than 

d) was chosen in the denominator as a scaling factor in 

max 

order to assure relatively good convergence also in areas of 
small (J)p values. For the average solution using the Two- 
Equation model of turbulence with a 41 x 24 grid system 
(Figure 4.4, page 97) this was normally achieved after about 
50 iterations, taking about 7-1/2 minutes on an IBM-360-65 
computer . 

IV. RESULTS AND DISCUSSION 


The results obtained -in the present study include a 
laminar flow solution and three turbulent flow solutions, 
the first of which uses the mixing length model, the second. 
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the turbulent energy hypothesis with a prescribed mixing 
length distribution, and the third, the Two-Equation model 
(all described in Section II of this chapter) . For the 
latter turbulence model a parametric study of two parameters 
of -the approach wind profile, uj and z^, was carried out. 
Before these results will be discussed, some light shall be 
shed on factors influencing the solution and affecting the 
accuracy from a different viewpoint than that discussed in 
the previous sections. It is referred to the empiricism 
involved in the solution procedure entering not only through 
a "proper" choice of coefficients and Cg^ in 

Equations 4.16, 4.19 and 4.23, but also through- the selec- 
tion of "suitable" boundary conditions whose dominating 
effects are not always recognized. 

Factors Influencing the Solution 

Boundary conditions . In many cases it is quite easy 
to conjure up some kind of plausible boundary conditions, 
but attempts to determine boundary conditions which are 
equally realistic, accurate and stable can be highly 
frustrating and often their selection ends with the compro- 
mise that the first condition is neglected in favor of the 
last two. 

This is even more the case when conditions cannot 
properly be formulated because the necessary empirical 
information is not available. Unfortunately, the only way 
around this problem is to extend the computations far enough 
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upstream or downstream that either realistic assumptions can 
be made (vanishing gradients) or the influence from the 
boundary conditions becomes negligibly small for the region 
of interest in the flow field. 

This problem arises not so much with the formu- 
lation of the more familiar boundary conditions for vor- 
ticity and stream function, which, with the exceptions 
discussed earlier, are believed to be rather unequivocal, 
partially because one can draw heavily from wind tunnel 
tests or flow visualization experiments. 

The problem appears, however, with the more unknown 
variables k* and I*, with the formulation of the outflow 
boundary condition being the most controversial. Should the 
turbulence kinetic energy, for example, be allowed to decay 
back to the original free stream value (Equation 4.52) or 
should it exceed this value in dependence on the higher 
local friction velocity (Equation 4,46)? Should it be 
constant in the z*-direction ox should its streamwise 
variation be zero? Fortunately, none of the above conditions 
posed any convergence problems so they could all be investi- 
gated. The results are shown in Figures 4.18 and 4.19. It 
should be mentioned that the stream function as well as the 
vorticity pattern were hardly affected by the above changes 
in the k*-boundary condition and are therefore omitted. A 
comparison between the turbulence kinetic energy distri- 
butions of the first two cases shows that the Influence of 
the outflow boundary conditions 
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Figure 4.18. Influence of outflow boundary condition on k-solution; 
(a) k* = (u*/Cyj)2 and (b) (3k*/9x*) = 0. 
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Figure 4.18 


(continued) 
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(a) 


k* “ 



(4.52) 


and (b) 


3k* 


0 


is limited to a far dovmstream region about seven step 
heights behind the step. Apart from this region the two 
solutions are the same and, therefore, independent of their 
particular boundary condition. The third case with 


k* 


II ★ 

"*out 



(4.46) 


shows a somewhat different picture. The changed brought 
about by this condition seem bigger and reach farther into 
the flow field than in the previous case. While with the 
first two conditions the k*-decay in z*-direction is faster 
before the step than it is behind it, the revfsrse is true 
here. This behavior is partially a result of the upper 
boundary condition 3k*/3z* * 0. More realistically this 
boundary condition should be 3k*/3x* = 0 which is equivalent 
to k* = const. Substituting this boundary condition into 
the case presented in Figure 4 .18 (a), one obtains the dis- 
tribution shown in Figure 4.20. Comparing the two one finds 
hardly any or little changes. 
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Figure 4.20. Influence of upper boundary condition on k-solution, 
Ok*/3x*) » 0. 
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In conclusion/ it can be said that In the present 
study, except for the conditions chosen in Figure 4.19, page 
133, the upper and outflow boundary conditions for k* do not 
significantly influence the solution. The question if (a) 
or (b) , i.e,, k* = const, or k* = f(z*), is the more 
realistic boundary condition becomes superfluous, because 
apart from the immediate vicinity of the boundary they gave 
the same results. For this reason and the fact that pre- 
scribed boundary’ values give generally better convergence 
than normal gradient type conditions, the k*-boundary con- 
ditions of Figure 4.20 were mainly used in the following 
computations. Similar arguments hold for the boundary 
conditions of the turbulence length scale A*., 

Empirical coefficients . It was already pointed out 
earlier that the empirical constants appearing in Equations 
4.19 and 4.23 are no universal constants, but depend on the 
particular flow case under consideration. It is common 
practice to determine these constants in a preliminary 
evaluation by applying the governing equations to simple 
flow situations such as Couette flows 135] , flows with 
homogeneous turbulence behind a grrd [36] or turbulent 
boundary layers in local equilibrium where the generation of 
turbulence at any point in the flow field is balanced by the 
local dissipation [36] . To achieve closure, the relations 
found this way will then have to be supplemented with infor- 
mation gained from experimental data. However, if 
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experimental results are not available, there is little hope 
for an exact determination of the constants and the above 
path has to be abandoned. 

Instead, we shall proceed in a different fashion and 
shall look at the individual constants and examine how their 
variation effects the solution. Also, a trial set of 
"universal" constants deduced from an early proposal of 
Spalding [43] for a Two Equation model actually using length 
scale as the second turbulent transport equation will be 
tested. Finally, a set of preliminary constants will be 
selected as a result of the foregoing investigation. 

Table 4c 2 gives a summary of the individual sets of 
constants reviewed. Case 1 is the reference solution. 

Cases 2 through 6 represent attempts to assess the 
respective influence of variations in the effective vis- 
cosity (ACyj,), in the turbulence kinetic energy dissipation 
rate (ACd^), in the effects of length scale stretching 
(ACsfl) and of the length scale breaking (AC 3 ^) on the 
solution. The investigation included the survey of the 
respective ip*, k* and distributions. 

The results are shown in Figures 4.21 through 4.28. 
The turbulence kinetic energy distributions have been non- 
dimensionalized with the turbulence kinetic energy of the 
undisturbed approach flow and for easier comparison the 
turbulence length scale variations for the different cases 
are presented in a single figure (4.28) where only the 
respective lines of constant passing through one and the 
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TABLE 4.2 

SET OF EMPIRICAL CONSTANTS USED IN TRIAL CALCULATIONS 


Case 

Cyo 

CDo 

CSo 

CBo 

1 

1.0 

1.0 

1.0 

1.0 

2 

0.1 

1.0 

1.0 

1.0 

3 

1.0 

0.1 

1.0 

1.0 

4 

1.0 

1.0 

0.1 

1.0 

5 

1.0 

1.0 

1.0 

0.1 

6 

2.0 

l.p 

1.0 

1.0 

7 

1.74 

1.18 

0.27 

0.14 

8 

1.0 

0.1 

0.27 

0.14 
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Figure 4 •21, Streamline pattern and turbulence kinetic 
energy contours for case 1 (Cy^ == 1.0, Cp. = 1.0, 
^So * 1*0, = 1.0. 
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Figure 4.23. Streamline pattern and turbulence kinetic 
energy contours for case 5 (C,.„ = 1 . 0 , Cn =1.0. 
Cs, =• 1.0, Cb, - 0.1). 
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Figure 4.27. Length scale distribution for case 1 (lines of constant St*). 
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Figure 4.28. Lines of constant length 
the different cases. 
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same upstream point P are plotted. Case 1^ the reference 
solution, with all constants set equal to unity shall not be 
further discussed at the moment. We shall return to this 
solution later. 

In case 2 with Cy, = 0.1 for which no figure is 
shown, it was intended in view of Equation 4.16 and the 
first term on the right side of Equation 4.19 to cut down 
the turbulence kinetic energy production. This was easily 
achieved, however, the amount of reduction was too large. 
Everywhere in the flow field k* was smaller than the 
approach condition of unity, implying that the step would 
reduce turbulence. A reasonable reduction should be 
obtained for a value of Cy^ only slightly smaller than one. 

An increase in k* was obtained in case 6 (Figure 
4.24) with Cyij “ 2.0. Compared with the reference solution, 
a doubling of Cy^ produced a maximum k* about four times as 
large. 

The turbulence kinetic energy levels should also 
increase if the dissipation rate is diminished, thereby 
increasing the net-production (Equation 4.19). This was 
done in case 3 by setting C^^ ® 0.1. The resulting k*- 
distribution (Figure 4.22) has a maximum of about 3-1/4 
times that of the reference case, which like case 6 is 
unrealistically high. 

The attempt to decrease the length scale stretching 
contribution (Equation 4.23) by reducing Cg^ to 0.1 in case 
4 failed, because the numerical procedure diverged. 
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However, decreasing the length scale breaking con- 
tribution with Cgg = 0.1, thus increasing the net H*- 
production and the k*-generation (case 5) had a similar 
effect to decreasing the dissipation in case 3. The k*- 
levels rose in the maximum to about 3-1/2 times the refer- 
ence values. 

Thus, there are four effective means to increase the 
turbulence kinetic energy levels; 

1. Increase 

2. Increase CSg. 

3o Decrease 

4. Decrease Cb^. 

Reducing the turbulence kinetic energy requires opposite 
measures but from the experience with cases 2 and 4 we know 
that these have to be applied with greater moderation. 

In the above trials only one of the respective 
constants was changed at a time to test the individual 
effect, but due to the coupling of the governing equations, 
changes of several parameters at a time can have quite 

m 

different effects than can be estimated from the individual 
behavior and a more systematic investigation is required. 

However, we shall restrict ourselves in the current 
study to the above cases and will only, out of curiosity, 
try Spaldings [43] coefficients (case 7). The resulting k*- 
distribution is as high as some of the previous cases 
(Cy^ = 2.0). This is probably a result of the relatively 
large value for Cy ^ . 
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Other researchers from the same group of Gosman, 
et al. [20] at the Imperial College [36# 37] later used much 
lower values for the constants in the k*-equation and these 
values were also tried, case 8, The constants for the l*- 
equat ion from [43], i,e., Csg = 0.27 and Cbj, = 0.08 were 
still retained, however. Figure 4.26, page 144, shows the 
corresponding k*-distribution, which surprisingly enough is 
almost identical to the one in case 3 where and had 

the same values,- though different coefficients for the i*- 
equation were used. It seems that the influence of Cy^ and 
CDq dominates the influence of the coefficients in the 
length scale equation. 

Until this point we have only looked at the turbu- 
lence kinetic energy distributions, but the goodness or 
physical reality of a solution is probably better recognized 
from the stream function distribution or even the length 
scale plot. 

It is interesting to note that in the respective 
stream function plots there are hardly remarkable changes, 
except for the cases 1 and 6 which are the only ones pro- 
ducing a noticeable separation from the corner. 

It is, therefore, not without reason that the length 
scale distributions in Figure 4.28, page 146, for these two 
cases are different from those of the remaining cases. They 
are not only different, but more important physically more 
meaningful. This insofar as it seems quite correct that %* 
decreases as one approaches the step parallel to the wall. 
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due to the acceleration in the flow field, stretching the 
vortices or due to the breaking of the larger eddies by the 
increasing shear. 

It can be seen in Figure 4.27, page 145, that closer 
to the wall the reverse is true. Here the increase of A* 
along a line of constant height represents the dissipation 
of energy from smaller eddies and the formation of larger 
eddies due to the deceleration of the flow which finally 
results in the formation of the front separation region. 

A similar phenomenon should occur near the upper 
separation region in those cases where separation occurs 
(cases 1 and 6) . It seems realistic that after the region 
of decreasing A* in front of the step a region of growth of 
A* should develop because it is essential for the top 
separation or recirculation region to form that energy is 
dissipated from the smaller eddies into this region of 
larger A*. 

This estimated correct behavior in the A*-distribution 
is predicted only in cases 1 and 6, which from this view- 
point seem to be the most realistic. They are also the only 
ones to create a top separation region in the t|;*-plot which 
is experimentally known. In case 6, however, the turbulence 
kinetic energy distribution was very high, and not as 
realistic as for the reference case. 

For these reasons it is believed that case 1, where 
all "empirical" constants are set equal to unity, presents a 
physically meaningful description of the flow problem under 



151 


investigation. Therefore, it is recommended that until 
suitable experimental results are available to redefine the 
values of the individual constants, reasonable realistic 
predictions can be expected with values of unity for these 
constants, which were used throughout all subsequent 
calculations . 

Comparison of Different Models 

Before looking at the results of the three turbu- 
lence models we shall take a brief look at some results of a 
"laminar" solution, which were obtained in the early stages 
of the progrcim checkout, when the computer code for turbu- 
lent flow with the appropriate boundary conditions was 
tested with the turbulent viscosity set equal to zero, to 
save computing time. 

Laminar solution . As already indicated above we are 
not dealing with a truly laminar solution because of two 
reasons: First, the inflow and outflow boundary conditions 

for the stream function and the vorticity were based on the 
s^une logarithmic approach velocity profile used in the 
subsequent turbulent cases. A parabolic power law profile 
would have been a more suitable condition. Second, the 
occurring Reynolds number based on the step height 

u*h* 

Re„ - (4.72) 


where 
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u! 


u* = 


Zn 


h* 


z*o 


( 4 . 73 ) 


of 3.93 10® lies certainly outside the laminar flow regime. 

Therefore the "laminar" solution would probably better be ^ 
called "quasi-laminar , " in the sense that only the laminar 
or molecular viscosity is used for the calculation of a flow 
field which should actually be turbulent. Although one may 
question the validity of the solution under these circxim- 
stances, it exhibits some interesting features and shows 
typical laminar behavior when compared to some of the later 
turbulent solutions. 

Figures 4.29 and 4.30 show the streamline distri- 

\ 

bution and a plot of the velocity profiles in the vicinity 
of the step. The laminar character of the flow can best be 
recognized by taking a short look at one of the subsequent 
turbulent solutions^ for example Figure 4.31, and comparing 
the two solutions. We know from experience that laminar 
flow favors separation, while turbulence tends to suppress 
it, or leads to fast reattachment of the flow. The same 
trend can be observed here. In the laminar solution the 
flow ahead of the step separates earlier than in the turbu- 
lent case. In general, the disturbance caused by the step, 
recognized by the deflection of the streamlines, is felt 
much further upstream by the laminar flow than it is by the 
turbulent flow. While separation from the upper step corner 
is strong in the laminar case with the separation region 
extending far downstream, reattachment occurs already a 
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.30. Velocity profiles for laminar flow solution. 
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Figure 4.31. Streamline pattern for turbulent flow solution using the 
Prandtl-mixing-length concept (Re^ = 3.9 x 10 ®). 
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short distance behind the upper step corner in the turbulent 
case. 

Comparison of different turbulent solutions . Results 
for the three different turbulence models employed are shown 
in Figures 4.31 through 4.36. A streamline pattern for the 
PML model, a streamline distribution and the turbulence 
kinetic energy contours for the TKE model and finally a 
streamline pattern with turbulence kinetic energy and length 
scale distributions for the Two-Equation model. 

From the respective streamline distributions it 
seems at first glance that all three solutions give approxi- 
mately the same answers. However, this is deceiving as 

% 

inspection of the remaining figures indicates and it is 
therefore appropriate to mention that the streamline pattern 
is not really a sensitive indicator of the correctness of 
the turbulence model. It is thus more informative to look 
at the turbulent quantities like turbulence kinetic energy, 
shear stress or effective viscosity, where the, last two can 
readily be compared in all three solutions. 

In Figures 4o37 and 4.38 we see the individual 
effective viscosity predictions along the upper step wall 
for the near wall nodes as well as for the surface itself. 
Although measurements of shear stress near reattachment are 
not completely reliable, it is fairly well known that the 
surface shear stress 
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Figure 4.32. Streamline pattern for turbulent flow solution using the 
turbulence kinetic energy concept with prescribed length scale; 
One-Equation model (Re = 3.9 x 10®). 
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Figure 4.35. Turbulence kinetic energy contours for Two-Equation model. 
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Figure 4.37. Effective viscosity variation along near wall 
nodes on the top of the step for the different 
turbulence models. 
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Figure 4.38. Effective viscosity variation along the wall 
on the top of the step for different turbulence models. 
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and thus with it the effective viscosity rises rapidly after 
reattachmento However, as near the reattachment point the 
mean velocity gradients are low, the PML model with 
Equation 4,13 likewise predicts low values for or u*. 

Surprisingly enough, the results of the TKE model are little 
better and only the Two-Equation model follows the expected 
trend. 

The weakness of the PML model with Equation 4.13 
lies not only in the incorrect implication that p* vanishes 
whenever the mean velocity gradients are zero, but more 
generally in the fact that y* is assumed to 'depend only on 
local flow properties. We know, however, from experience 
that the local level of velocity fluctuations is determined 
not only by the events of the point in question, but also by 
influences which have originated some distance upstream. 

For instance, the high level of velocity fluctuations at the 
above reattachment point originates from a highly turbulent 
shear layer issuing from the upper step corner. It is 
important that this nonlocal character of turbulence is 
taken into account. 

The fact that the TKE model incorporates an addi- 
tional transport equation for k* permits account to be taken 
of the influence of neighboring regions on the local turbu- 
lence energy. However, as seen in Figure 4.37 and also 
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reported by [44] , a TKE model without a length scale trans- 
port equation produces results that for wall boundary layers 
are only insignificantly better than those obtained with 
PML models. To show how difficult it is to make an 
accurate guess for the length scale distribution entering 
the TKE model r reference is made to an Jl*-distribution 
actually calculated with the Two-Equation model in Figure 
4.36, page 161. It is the convective transport of Jl* which 
plays an important role in the determination of the flow 
properties and which results in more physically correct pre- 
dictions with the Two-Equation model. 

Taking a second look at the respective stream 
function plots one notices a difference in the geometry of 
the upstream and downstream separation bubbles for the 
different models. The predicted dimensions decrease from 
PML to TKE to Two-Equation model. Table 4.3 gives the 
approximate values for the individual cases. This decrease 
certainly is associated with the higher turbulence levels of 
the respective models. 

Generally the top separation region seems somewhat 
smaller in its vertical extent than one would expect. 

Besides this possibly being the result of the problem men- 
tioned in Section III of this chapter, it also seems logical 
that in an atmospheric boundary layer where the ratio of the 
step height to boundary layer thickness, i.e., h*/6* is 
small, smaller separation regions should occur than those 
reported for most wind tunnel experiments where h*/6* is 
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TABLE 4.3 


GEOMETRIES 

OP SEPARATION REGIONS 
TURBULENCE MODELS 

FOR DIFFERENT 


MO(kel 

*S 

yR 

*R 

PML 

1.2 

0.6 

2.2 

TXE 

1.0 

0.6 

1.4 

Two-Equation 

0,6 

0.4 

1.2 




167 


usually large. More on this topic will be found in the 
parametric study given in the next section. 

Let us now turn to a closer examination of the tur- 
bulence kinetic energy contours for the TKE and Two-Eguatibn 
model in Figures 4.33 and 4.35, pages 158 and 160. These 
can readily be related to the more commonly used turbulence 
intensities 




(4.74) 


by 



(4.75) 


In the Two-Equation model the energy contours reveal quite 
realistically that the shear layer, which grows from the 
corner of the step, generates high energy levels. These are 
carried downstream giving rise to comparatively high turbu- 
lence intensities in the downstream region of the reattach- 
ment point. This downstream convection of energy together 
with the increase in the level of length scale (Figure 4.36, 
page 161) conspire to produce particularly high turbulent 
viscosity in the region behind reattachment. 

This phenomenon is also visible in the TKE model k*- 
contours. However, the intensities in the downstreeim 
reattachment region are considerably smaller because of the 
inadequate prescription of the length scale, thus resulting 
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in the low y* predictions there, as already seen in Figure 
4.37, page 162. 

Another yet more eye-catching difference is the dif- 
fusion of k* in the TKE model as compared to the Two- 
Equation model, especially in the downstream flow region. 
Experience tells us that turbulent shear layers usually 
spread at a rate x^ with m depending on the type of flow. 
Likewise should the shear layer originating from the step 
corner spread out and with it the turbulence kinetic energy 
or turbulence intensity. Again this behavior is only in the 
Two-Equation model described correctly. 

In the flow region ahead of the step the turbulent 
intensities do not differ that greatly for the two models. 
Qualitatively the predictions seem to agree quite well with 
the experimental findings of Taulbee and Robertson [45, 46] 
who report that the turbulent intensities increase strongly 
toward the front separation point. Primarily there is a 
spreading of the zone over which the higher intensities 
extend rather than a significant increase in magnitude. The 
peak intensities occur appreciably outside the dividing 
streamline, which defines the outer edge of the front sepa- 
ration bubble. Outside this peak the turbulence decreases 
with z* as it does in any boundary layer. Inside the sepa- 
ration bubble the turbulence decreases with distance to an 
essentially small constant value. The authors point out, 
however, that the peak intensities and shear levels reached 
here are still considerably less than values typical of 
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those in fully developed jets or free shear layers, like the 
one originating from the step corner (compare with Figure 
4.35, page 160) . 

Finally, to complete the presentation of the 
dependent variables involved in the calculation, a vorticity 
contour plot determined from the Two-Equation model is given 
in Figure 4.39. It shows very clearly the downstream 
influence of the step on the vorticity pattern, the creation 
of a highly disturbed region with steep vorticity gradients 
not only normally to the wall but also in the flow direction. 
The outer edge of this disturbance or wake region can be 
approximated as indicated by the dotted line, which, for 
this particular flow case, is roughly of the form z* ® 
(x*/h*)°‘®. 

It is also interesting to note that apart from this 
region and the immediate upstream vicinity of the step, the 
constant vorticity lines follow quite closely the streamline 
pattern so that the application of the "Frozen Vorticity 
Theory” of Taulbee and Robertson (45, 46] to the current 
flow problem seems to be a reasonable assumption in this 
flow region. "Frozen Vorticity" assumes that for certain 
flow situations where boundary changes and effects occur so 
suddenly that viscosity or turbulent mixing do not have any 
time to act on the flow, the diffusion of vorticity can be 
neglected. The vorticity is then fluid bound, i.e=, 

"frozen" to the stream function and its distribution fixed 
as given by the approach flow. The governing stream 
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Figure 4.39. Vorticity distribution for Two-Equation model 
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function equation (4.10) would reduce under these circum- 
stances to 

(4.76) 

To summarize the )cey results of this paragraph, one 
might say that all three turbulence models, PML, TKE, and the 
Two-Equation model, yield quite realistic streamline patterns 
and velocity distributions. In the prediction of turbulent 
flow properties like k* and t* the Two-Equation model was 
found to be superior. The inadequacies of the PML and TKE 
model are that they either neglect (PML) or do not ade- 
quately account for (TKE) the convective transport of 
turbulence quantities which in the current flow problem is 
especially important as the step strongly influences the 
downstream flow region. 

Effect of Variation of the Parameters of the 
Approaching Wind 

A parametric study was conducted utilizing the Two- 
Equation model to assess the influence of the parameters of 
the approaching wind profile, u$ and z*, as well as the step 
height h* on the solution. 

Before proceeding to the calculation and discussion 
of selected flow solutions it is useful to study the 
dimensionless aspects of the flow problem. 

A nondimensional form of tne governing equations was 
obtained by adopting the step height h* as the characteristic 
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length and the friction velocity u* as the characteristic 
velocity. 

Introducing 


0 ) 


(o*h* 



uj ' 

p*ujh* 

k* 

0 

z* 

uJ2 '■ 

X, 

h* 

- . 


z* 

h* ' 

z 

h* 


- 2 * 

2o - ij* ; 


Re* = 


p*uth* 


u* 


( 4 . 77 ) 


The effective viscosity can then be nondimension- 
alized as 


eff p*u*h* Re. 


V] 


p*ufh* 


( 4 . 78 ) 


which yields for the Prandtl-Kolmogorov formulation of the 
turbulent viscosity 


'eff 


1 

Re* 


+ 


1 / 2 , 


( 4 . 79 ) 


The resulting nondimensional form of the governing equations 
then is: 

Stream function equation (continuity) 



( 4 . 80 ) 
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Vorticity transport equation (momentum) 
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Turbulence kinetic energy transport 
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and the transport equation for the turbulence length scale 
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(4.83a) 


with 




(4.83b) 


After nondimens ionali zing also the boundary conditions, 
e.g., the inlet conditions in the four variables rl>, o), k and 


Z: 
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^ = i| 

K ' 


(z + 2 o) |jln — - ij + zo 


(I) = - 


1 

K (Z + Zo ) 9x^ 


k = 1 

I = K (Z + Zo ) 


(4.84) 


one finds that the problem reduces from the three 
dimensional parameters z* and h* to the two dimensionless 
groups Zo and Re^. 

However, in view of Equation 4.72 and the high 
Reynolds number usually occurring in atmospheric flows, 
where 



(4.85) 


the flow problem becomes essentially independent of Reynolds 
number. Thus, the only remaining significant parameter is z©. 

In the subsequent study calculations were carried 
out for the following five values of Zo. 

(1) Zo = 0.005 

(2) Zo = 0.020 

(3) Zo = 0.045 

(4) Zo = 0.075 

(5) Zo = OclOO 

The resulting k- and )l-distributions for the 
cases 1, 3 and 5 are shown in Figures 4.40 through 4.42. 
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.41. Turbulence kinetic energy contours for (a) zo = 
Zo = 0.045, (c) Zo = 0.1, 
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Figure 4.41. (continued) 
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Comparing the various streamline patterns one 
notices that the shapes of the upstream and downstream 
separation regions vary with Zq. While the top reattachment 
length changes substantially, the front separation 
distance Xg varies only slightly. Figures 4.43 and 4.44 
reveal the geometry of the separations and reattachments as 
characterized by the distances x_, y_ and x_. 

The reattachment length was found to increase 
with decreasing Zo, indicating that a smoother surface or 
larger step height would delay reattachment of the flow on 
top of the step. The reason for this is readily understood 
by looking at the respective TKE- and 5, -contour plots in 
Figures 4.41 and 4.42. For smaller Zo the 'turbulence 
intensities and length scale are seen to decrease in the up- 
stream and downstream vicinity of the upper step corner. 

This is a result of the higher flow acceleration created by 
the displacement of a fuller approach velocity profile con- 
taining higher momentum near the wall (see also Figure 
4.45(b)). Acceleration in the flow is generally known to 
diminish turbulence production and to reduce the turbulence 
length scale or typical eddy size by means of vortex 
stretching. These lower turbulence intensities together 
with the reduced length scale lead to a lower effective 
viscosity or shear on top of the step, delaying reattachment 
and causing Xj^ to increase. 

For the forward separation region the situation is 
somewhat different. Except for a small region. in the lower 
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Figure 4.44. Forward separation region; (a) geometry in 
dependence of zo, (b) maximum vorticity in dependence 
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Figure 4.45. Velocity profiles for various Zo; (a) approach 
velocities, (b) velocity overshoot at x = 0 due to the 
step for different zo. 
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Zo range ^ the bubble is found to grow in size with in- 
creasing surface roughness parameter, even though the 
maximum vorticity occurring inside the bubble decreases, as 
shown in Figure 4.44(b). In the zo -range under conside- 
ration where the vorticity decreases only slightly the 
growth of the bubble seems to be influenced largely by the 
growth of the turbulence length scale. In the lower zo- 
range mentioned above, however, the steep decline of 
vorticity seems to be more influencial than the growth of 
the length scale resulting in a partial reduction of bubble 
size. 

Taulbee and Robertson [46] presented the results of 
their investigation for smooth walls in dependence of the 
ratio of step height to boundary layer thickness h*/6*. 

This ratio is small for most atmospheric boundary layer 
flows and changes only little; however, it is similar to the 
ratio h*/z* used here, in the sense that both parameters 6 * 
and z* characterize the nature of the approaching velocity 
profile in relation to the height of the obstruction. 
Although a direct comparison with the above results is not 
possible because the functional relation between 6* and z* 
is unknown, one can nevertheless compare common trends. 
Taulbee and Robertson's results predict for 

^ = 1.8 = 3.8 


and for a step height a fourth of the above; 
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h*/4 _ 


0.45 Xj^ 


2.3 


For the present study the surface roughness of Zo = 0.005 
predicts the same reattachment length as h*/6* = 1.8, i.e., 

* 

gf = 0.005 -*■ Xj^ = 3.8 

and for a fourth of the step height we get 
= 0.02 = 2.35 

Carrying on the comparison in this fashion for several other 
step heights, one finds the predictions for the reattachment 
length Xj^ in almost exact agreement with the results of [46] . 
Proceeding the same way with dimensions Xg and yj^ of the 
front separation bubble, however, agreement is not as good, 
with the predictions of the Two-Equation model being con- 
sistently lower than Taulbee and Robertson's results. As 
already noted during the comparison of the three turbulence 
models (see Table 4.2, page 138) this is possibly a con- 
sequence of the higher turbulence and turbulent viscosity in 
the k-H-model and could be adjusted by a redetermination of 
the empirical constants. 

More information is contained in the gradients of 
the stream function, i.e., the velocity distributions. 

Figure 4.45, page 188, indicates how the presence of the 
step affects the horizontal velocity profile at x * 0.0 for 
the different roughness parameters. The profile with the 
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strongest velocity gradient at the wall and the highest 
momentum near the surface produces the largest velocity 
overshoot. This is easily perceived because for the smaller 
Zo the flow carries more mass near the wall, which is 
suddenly displaced by the step, creating a locally acceler^ 
ated flow region. This result was also obtained in the 
boundary layer analysis. The overshoot is in accordance 
with the lower shear predictions and longer reattachment 
lengths for the smaller z g . 

For completeness. Figure 4.46 shows a set of 
velocity profiles over the entire flow region for case 1 
with Zo = 0.005. 

The turbulence kinetic energy contours reveal that 
for increasing zo the turbulence energy levels in the shear 
layer growing from the step corner rise accordingly. Like- 
wise, the turbulent viscosity and with it the turbulent 
shear stress increases in unison. The same trend becomes 
apparent by looking at Figure 4.47 which, for easier 
comparison, shows the individual turbulence profiles 
together at the x «• 3.0 downstream location. Once more this 
gives the reason, why for smaller surface roughnesses bigger 
overshoot and larger downstream separation regions are 
expected. 

The following Figure 4.48 gives again a complete set 
of turbulence intensity profiles over the whole flow field 
for case 1 with Zo - 0.005. As already pointed out in the 
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previous sections, these profiles essentially agree with the 

findings of Taulbee and Robertson [46] . | 

I 

The turbulence length scale plots expose; how 
different the ii- -distributions actually can become frbm the 
linear variation originally assumed in Equations 4.36 and 
4»47. The downstream influence of the step on the turbii- 
lence stf'udture iS recognized by the growth in SL indicating 
an increase of typical eddy size in the spreading Shear 
layer o To facilitate a comparison between the various 
roughness parameters and their effects on the turbulence 
structure. Figure 4,49 shows the -profiles at x = 3.0 for 
the different Zg. The'profiles have been nondimenSiOnalized 
with Equation 4. 47 to allow easy comparison with ; the linear 

distributiono In accordance with earlier results thS 

j 

rougher surfaces produce shear layers with sizeable in- 
creases in the typical eddy size* which through Equation 
4.16 au^ent' the turbulent viscosity and tend to ^shorten the 
upper serration .regions o Pot small zo , however,- the H pre- 
dictions, are below those given in Equation 4. 47 hence 
leading to higher overshoot and delayed reattachment. 

Prom ’the four variables oo, tp, k and the vorticity 
and the zo -influence on it has not yet been discussed. 
Recall, Figure 4.39, page 170, showed th^ vorticity contours 
for zo - 0o045. It was pointed out that a disturbance 
region with locally high streamwise vortrcity gradients is 
formed and spreads downstream parabolically. Figure 4.50 
shows how a surface roughness change affects the extent of 
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this region of high vorticity. It is seen that a Zo- 
increase generates a wider spread, possibly as a consequence 
of the changes in the turbulence structure (larger eddies) 
of the shear layer, discussed in the previous paragraph. 

Concluding this parametric investigation, one may 
summarize the results by stating that in the downstream flow 
region an increase in Zo gives rise to higher turbulence 
levels in the shear layer originating from the step corner. 
This, in turn, results in higher shear stress leading to 
fast reattachment. The typical eddy size in the shear layer 
increases with zo and causes rapid spreading of a region of 
high streamwise vorticity gradients. 

In the upstream flow region changes in Zo have only 
moderate influence on the flow parameters. Except for very 
small surface roughnesses, the front separation bubble grows 
in size for increasing z^. 
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